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A  COMPARISON  OF  SOME  SIGNAL-PROCESSING  ALGORITHMS 
TO  SUPPRESS  TOW-VESSEL  NOISE  IN  A  TOWED  ARRAY,  WITH  RESULTS  FROM  A 

SHALLOW-WATER  SEA  TRIAL 


by 

Ulrich  E.  Rupe  and  Even  B.  Lunde 


ABSTRACT 


A  short  trial  has  been  conducted  in  shallow  water  using  a  towed  array  and  a 
number  of  signal -processing  techniques.  Results  are  presented  on  the 
outputs  from  conventional  beamforming  and  maximum-likelihood  beamforming 
and  on  some  basic  techniques  for  suppressing  towship  noise.  The  signal¬ 
processing  package  is  implemented  in  software  and  provides  the  various 
forms  of  processing  by  using  the  array  covariance  matrix  and  eigenvalue/ 
vector  analysis.  Possible  future  work  on  these  topics  is  proposed. 


INTRODUCTION 


The  self-noise  of  the  platform  can  enter  a  hull -mounted  sonar  array  from 
different  directions  and  in  various  frequency  bands,  thereby  presenting  a 
number  of  problems  in  the  reception  of  signals  during  the  initial  processes 
of  target  detection,  classification,  and  tracking  when  signal  levels  are 
low.  The  advent  of  the  towed  array  has  minimized  these  problems 

considerably  and,  in  essence,  reduced  them  to  simply  one  of  a  noise  source 
ahead  of  the  array. 

However,  depending  on  the  characteristics  of  the  towed-array  system  (the 
scope  of  the  tow  cable,  the  beamwidth  of  the  array,  etc.)  and  the  environ¬ 
mental  conditions  in  the  area  of  operation  (depth  of  water,  type  of  sea 
bottom,  etc.),  the  towship  noise  can  appear  in  a  number  of  beams  in  the 
ahead  ±90°  sector  and  cause  serious  operational  limitations. 

This  memorandum  compares  a  number  of  signal-processing  techniques  by 
describing  results  obtained  from  a  sea-trial  in  shallow  water  using  a 

towed-array  system. 

The  general  headings  are: 

a.  The  presence  of  a  number  of  signal  paths  from  the  towship  to  the 
array. 

b.  The  effect  of  the  basic  techniques  for  suppressing  towship  noise. 

c.  Target  indication  and  bearing  estimation. 
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FIG.  1  SEA  TRIALS 

a)  Runs  1  and  2 

b)  Run  3 
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FIG.  1  SEA  TRIALS 

c)  Run  4 

d)  Run  5 
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The  conclusions  and  recommendations  outline  the  proposals  and  arrangements 
for  further  studies  and  sea-trials  with  the  aim  of  reducing  the  operational 
limitations  caused  by  towship  noise  and  improving  the  estimation  of  target 
bearing. 


1  METHOD 
1. 1  Sea  Trials 

Sea-trials  were  conducted  in  the  area  off  the  northwestern  Italian  coast 
indicated  on  Fig.  1.  Some  parts  of  the  area  have  a  constant  water  depth  of 
150  m;  in  other  parts  the  depth  increases  from  150  m  to  as  much  as  500  m. 
Bottom  data  are  assumed  to  be  the  same  as  those  measured  by  SACLANTCEN  in 
an  area  a  few  kilometres  further  south  (Area  D1  of  <1>).  A  Zh  m  core 
(Core  127)  in  that  area  was  mainly  of  silt  with  some  shells,  overlain  in 
the  upper  50  cm  by  sand  and  then  by  clay.  The  relative  sound  speed  was 
1.025  and  the  relative  density  was  about  1.8.  No  data  on  bottom-reflection 
losses  are  available.  For  the  theoretical  predictions  of  transmission  loss 
in  the  experimental  area,  sound-speed  data  have  been  taken  from  <2>. 

During  the  trials  the  sea-state  was  SS-0  most  of  the  time  and  never  higher 
than  SS-1.  The  towed  array  was  deployed  from  SACLANTCEN1 s  research  vessel 
MARIA  PAOLINA  G.  (MPG)  as  shown  in  Fig.  2  and  most  of  the  sea-trials  were 
conducted  at  a  speed  of  5  kn.  SACLANTCEN' s  workboat  MANNING  transmitted 
334  Hz  signals  from  a  source  suspended  at  depths  of  either  50  or  100  m. 
These  signals  were  used  primarily  as  a  CW  source  with  different  source 
1 evel s. 

Figure  1  shows  the  relative  positions  of  the  MARIA  PAOLINA  G.  and  MANNING 
for  the  five  runs  constituting  the  trials.  In  Runs  1  to  4  (Figs,  la,  b,  c) 
the  MANNING  was  at  anchor  whilst  transmitting;  during  Run  5  (Fig.  Id)  it 
was  underway.  The  various  runs  placed  the  MANNING  on  all  possible  bearings 
and  at  ranges  from  1  km  to  30  km.  Figure  3  gives  the  sound- speed  profiles 
recorded  during  the  sea- trial.  They  show  that  a  shallow  sound  channel 
usually  existed  at  about  100  m. 


1. 2  Towed- Array  System 

SACLANTCEN ‘s  Prakla-Seismic  towed  array  was  used,  with  a  maximum  tow-cable 
scope  of  700  m.  The  hydrophone  spacing  employed  was  1.96  m,  corresponding 
to  A./2  at  approximately  385  Hz. 

An  additional  hydrophone,  referred  to  as  the  reference  hydrophone,  was  used 
to  investigate  one  of  the  basic  methods  of  suppressing  towship  noise.  This 
could  be  clamped  onto  the  tow-cable  at  distances  of  up  to  200  m  behind  the 
towship  (Fig.  2).  No  streamlined  body  was  available  for  this  additional 
hydrophone  during  this  sea  trial  and  heavy  flow  noise  was  expected  to 
degrade  the  correlation  with  the  array  signals.  The  signal  from  the 
reference  hydrophone  was  brought  in-board  for  use  in  the  signal -processing 
techniques  (see  Ch.  3). 
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FIG .  2  DEPLOYMENT  OF  TOWED  ARRAY  DURING  TRIALS 
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FIG.  3  SOUND-SPEED  PROFILES  DURING  TRIALS 
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1. 3  The  On-Board  Processing  and  Analysis  System 

A  block  diagram  of  the  on-board  processing  and  analysis  system  used  in  the 
sea-trials  is  shown  in  Fig.  4.  It  consisted  of  two  main  parts:  the  on¬ 
line  monitoring  system  and  the  off-line  processing  system,  as  described  in 
more  detail  in  the  following  chapters. 


2  ON-LINE  DATA  MONITORING 

The  on-line  data  monitoring  system  on  board  the  MPG  provided  a  display  of 
either  the  output  of  a  conventional  beamforming  system  (CBF)  or  of  the 
maximum- 1 i kel i hood  method  (MLM)  of  beamforming  for  a  single  selected 
frequency.  The  display,  which  used  the  standard  computer  line-printer  in  a 
novel  manner,  gave  an  immediate  visual  indication  of  the  power  in  the  CBF 
or  MLM  beams.  The  component  parts  of  the  system  are  described  below. 


2. 1  Signal  Conditioning  Unit 

The  lowest  cutoff  frequency  of  the  available  anti-aliasing  lowpass  filter 
bank  was  1  kHz.  As  the  frequency  range  of  interest  was  below  350  Hz,  a 
sampling  frequency  of  i  2.9  kHz  was  needed  to  have  an  aliasing  level  of 
=  -40  dB.  For  reasons  given  below,  the  chosen  sampling  frequency  of  the 
A/D  converter  was  12.288  kHz. 


2. 2  Digital  Lowpass  Filter  Bank 

The  Plessey  beamformer  was  used  as  a  lowpass  filter  bank  in  order  to  reduce 
the  sampling  frequency  by  implementing  48-point  linear-phase  finite- 
impulse-response  (FIR)  filters  with  3.072  kHz  (=  12.288/4  kHz)  input 
sampling  frequency,  and  819.2  Hz  (=  12.288/15  kHz)  output  sampling 
frequency.  The  beamformer  was  also  used  to  delay  the  data  from  the 
reference  hydrophone. 


2. 3  On-Line  Processing  and  Monitoring 

A  MAP  array  processor  and  an  HP  minicomputer  were  used  for  on-line 
processing  and  monitoring.  The  first  operation  was  a  1024-point  discrete 
fourier  transform  (DFT) ,  implemented  by  a  fast  fourier  transform  (FFT) , 
which  gave  a  resolution  about  0.8  Hz  (=  819.2/1024  Hz),  and  a  data  window 
of  1.25  s  length.  The  data  collection  and  FFT  processing  had  to  be  done  in 
sequence,  thereby  giving  an  "acquisition"  time  of  approximately  2.4  s.  In 
other  words,  a  time  window  of  1.25  s  was  repeated  each  2.4  s,  yielding  a 
time  gap  of  about  1.15  s.  The  DFT  outputs  for  ten  selected  frequencies 
were  written  on  magnetic  tape  for  use  by  the  off-line  signal -processing 
system. 

Flexibility  in  the  choice  of  the  parameters  for  the  calculation  of  the 
covariance  matrix  was  obtained  by  the  scheme  shown  in  Table  1.  The 
covariance  matrices  are  calculated  for  a  selected  number  of  frequencies 
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TABLE  1 

CALCULATION  OF  THE  NARROW-BAND  SPATIAL  COVARIANCE  MATRIX 


-  DO 

50  NS 

No.  of  averages 

-  DO 

20  NE 

No.  of  hydrophones 

CALL 

X(T) 

Time-domain  data  for  the 

phone  (NE)th 
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r-DO 
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from  filtered  and  fourier-transformed  time-domain  data  of  each  hydrophone 
channel.  The  required  average  time  is  obtained  by  repeating  this  operation. 


Thus  for  one  selected  frequency  out  of  the  ten,  the  output  vector: 


^1,S  ^2,S 


T 


was  used  to  build  up  an  estimate  of  the  covariance  matrix: 


R-  =  n  2  X(Mj+i)-X*(Mj+i)  ,  (Eq.  1) 

J  i=l 

where  N  is  the  number  of  hydrophones  used  to  build  up  the  output  vector 
X$  Tor  the  time  index  S  and  j  is  the  number  of  matrices  calculated. 

The  matrix  R.  was  used  to  produce  the  mean  power  outputs  for  181  beams, 

one  for  every  degree,  from  one  endfire  to  the  other.  These  outputs  were 
converted  to  decibel  levels  and  then  displayed  on  a  line  printer.  An 
example  for  conventional  beamforming  is  shown  in  Fig.  5.  The  bearing  is 
marked  above  the  top-line  for  a  narrow-band  processing  of  334  Hz.  All 
other  outputs  were  produced  with  the  off-line  processing  system.  Dependent 
on  the  number  of  acquired  time-series,  such  an  output  was  produced  for 
monitoring  purpose  at  intervals  of  1  to  2  minutes.  The  beamforming  was 
either  conventional: 


P  .  =  BkRA  k  =  1»  2»  •••181  (Eq.  2) 


or  "maximum-likelihood": 


k  =  1,  2,  ...181  ,  (Eq.  3) 


where  Bk  is  the  steering  vector  determining  the  beam  direction: 


,  exp  ( -  j  2tx  cos  pk  )  ,  . 1  T  ,  (Eq.  4) 


in  which  xn  is  the  coordinate  of  hydrophone  n  and  pk  is  beam  direction. 

Every  time  the  covariance  matrix  were  calculated,  the  "acquisition"  time 
was  increased  by  the  order  of  2  to  3  seconds.  Depending  on  the  number  of 
acquisitions,  the  beam  outputs  were  calculated  and  printed  from  the 
averaged  matrix.  The  minimum  number  of  averages  had  to  equal  the  number  of 
hydrophones. 
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CONVENTIONAL  BEAM-FORMING 

FREQUENCY  :  334  40 


e  0e  38.00  ae  00  pe  00  120  00  im  00  iso 


FIG.  5  EXAMPLE  OF  LINE-PRINTER  OUTPUT  FOR  CONVENTIONAL  BEAMFORMING 
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3  OFF-LINE  PROCESSING  AND  ANALYSIS  SYSTEM 

The  off-line  processing  and  analysis  system  used  tapes  of  discrete  fourier 
transform  (DFT)  data  produced  by  the  on-line  system  for  ten  preselected 
frequencies  and  gave,  as  outputs,  the  line-printer  display  mentioned  above 
for  the  on-line  monitoring  system  and  tabulated  data  relating  to  the 
various  parameters  of  interest.  The  off-line  analysis  provided  the  above 
outputs  for  the  following  signal-processing  techniques: 

a.  Conventional  beamforming  with  no  shading  (CBF). 

b.  Maximum- likelihood  method  of  beamforming  (MLB). 


c.  CBF  and  MLB  with  basic  techniques  for  suppressing  towship  noise 
by: 


1)  Use  of  a  reference  hydrophone. 

2)  Use  of  a  reference  beam. 

3)  Use  of  modified  eigenvalues. 


3. 1  General 

The  off-line  processing  system  consisted  of  an  HP  minicomputer  with  vector 
instruction  set  (VIS)  option  <3>,  a  magnetic  tape  station,  a  line  printer, 
and  a  terminal.  The  input  data  from  magnetic  tape  were  from  19  array 
hydrophones  and  one  reference  hydrophone,  quadrature  bandpass-filtered 
around  ten  selected  frequencies  and  with  a  bandwidth  of  0.8  Hz.  These  data 
were  then  processed  by  a  signal-processing  program,  and  the  output  written 
on  disc.  A  program  produced  plots  of  these  data  on  the  line  printer  in  the 
same  format  as  the  on-line  monitoring  program  (Fig.  5). 

In  addition  a  tabular  output  was  produced,  as  shown  in  Fig.  6,  for  each 

line-printer  output.  This  tabular  output  shows  a  whole  range  of  data;  for 

example,  columns  4,  13  and  22  give  the  beam  numbers  from  1  to  180,  columns 

5,  14  and  23  give  the  beam  angles  for  the  frequency  in  use,  and  the  columns 

under  CON  and  MLM  give  the  values  of  P.  .  (Eqs.  2  and  3)  for  the 

k  >  J 

conventional  and  "maximum-likelihood"  beam  processing. 

The  first  step  in  the  processing  program  is  to  create  the  following 
quantities: 


Ri=S  ^  X(Mj+i,q)-X*(Mj+i,q) 
j.q  ii  i=1 


(Eq.  5) 


(Eq.  6) 


ri  a  =  r  1  y(MJ+i«q)-y*(Mj+i>q) 

j,h  i=1 


(Eq.  7) 
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FIG.  6 


TABULAR  OUTPUT  CORRESPONDING  TO  LINE-PRINTER  OUTPUT 
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where 

X  :  Data  vector  from  array  hydrophones 

y  :  Data  scalar  from  reference  hydrophone 

R  :  Covariance  matrix 

Q  :  "Reference  vector" 

r  :  "Reference  scalar" 

q  :  Frequency  index 

i  :  Acquisition  time  index 

j  .  Time  index  after  averaging 

M  :  Number  of  acquisitions 


The  indices  are  omitted  hereafter. 

The  next  step  was  to  obtain  the  eigenvalues  and  eigenvectors  of  R.  As 
shown  in  App.  A,  R  can  be  written 


R  = 


N 

1 

n=l 


A 


n 


E  * 
n 


(Eq.  8) 


where  An  are  eigenvalues 
En  are  eigenvectors 
n  =  1.  2  N 

and  N  is  the  number  of  the  array  hydrophones 


There  were  three  reasons  for  introducing  the  "eigen-quantities".  Firstly, 
these  quantities  might  be  of  interest  on  their  own.  Several  signal¬ 
processing  methods  are  based  on  explicit  knowledge  of  these  quantities. 
Only  one  was  implemented  in  this  version  of  the  program,  the  so-called 
modified  eigenvalue  method. 

Secondly,  by  calculating  these  eigen-quantities,  one  might  gain  some 
feeling  for  their  behaviour  under  more  or  less  complex  sound-field 
conditions. 

Thirdly,  as  shown  in  App.  A  and  B,  the  ten  signal -processing  modes 
implemented  for  this  sea  trial  could  all  be  formulated  and  calculated  in 
terms  of  the  eigen-quantities.  The  signal  processing  in  the  off-line 
program  was  therefore  based  on  the  use  of  the  eigen-quantities,  whereas  the 
on-line  program  used  a  matrix  inversion  for  the  maximum- 1 ikel i hood 
technique.  A  check  on  the  two  programs  could  thus  be  made  by  comparing  the 
MLM  columns  of  the  outputs  illustrated  in  Fig.  6.  The  following  sections 
describe  these  methods  both  with  and  without  the  use  of  eigenvalues  and 
eigenvectors. 
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3. 2  Conventional  Beamforminq  (CB) 

The  mean  power  output  from  an  unweighted  (unshaded)  conventional  beam  can 
be  written  as  (Sect.  B.l  of  App.  B): 

P  =  B*R  B  . 

Introducing  eigenvalues  and  eigenvectors: 

N 

P  =  1  A  \  B*E  |2  .  (Eq.  9) 

n=l 


3. 3  "Maximum- Li kel i hood"  Beamforming  (MLB) 

The  mean  power  output  can  in  this  case  be  written  as  (Sect.  B.2  of  App.  B): 

p  =  (b^r-^)”1  ,  (Eq‘  10) 

and  with  eigenvalues  and  eigenvectors: 


(Eq.  11) 


3.4  Conventional  beamforming  plus  modified  eigenvalue  noise  suppression 

(CB  +  ME) 


The  eigenvalues  of  R  were  modified  according  to  some  rules.  One  such  rule 
was  to  make  eigenvalue  A.,  equal  to  the  smallest  eigenvalue  if  the 
corresponding  eigenvalue  E.  correlated  well  with  the  steering  vector 

representing  the  direction  of  the  towship.  Then: 

N 

P  =  Z  V  |B*Ej2  •  (Eq.  12) 


3.5  Maximum  likelihood  plus  modified  eigenvalue  (ML  +  ME) 

This  is  a  combination  of  ML  beamforming  with  modified  eigenvalues 

(Eq.  13) 


/  N  .n 
P  =|  I  A'  B*E 

l-i  1  n 
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3. 6  Conventional  beamforming  plus  reference  beam  noise  suppression 

(CB+RB) 

A  reference  beam  direction  with  steering  vector  B^  is  chosen  according  to 
some  rule  and  this  is  used  to  calculate  (Sect.  B.3  of  App.  B): 


=  (B*RB)  jjL- 1  B^RBr  |2  (B*RB  •  B*  R  Br)  J 


And  with  eigen-quantities: 


(Eq.  14) 


B*  R  B  =  I  X  I B^E  I 
n=l  n  1  nl 


(Eq.  15) 


N 

B*  R  B  =  Z  X  |B*E  I2 
r  r  _  n  I  r  n 
n=l 


(Eq.  16) 


N 

B*  R  B  =  I  X  (B*E)  (B*  E  )* 
r  ,  n  rn 

n=l 


(Eq.  17) 


3. 7  "Maximum  likelihood"  plus  reference  beam  noise  suppression  (ML  +  RB) 
Also,  from  Sect.  B.3  of  App.  B,  we  have: 


P  =  (B*R  1B)”1  |l-  |b*R  JBr 


12  _  1 

(B*R  B  •  B*R 
I  r 


V‘] 


(Eq.  18) 


And  with  eigen-quantities: 


B*R_1B 


I  X 
n=l 


_i 


B*E  !2 
n 


(Eq.  19) 


*  _  l 

B  R  B  = 
r  r 


* 

B  R  B 


Z  X  |B  *E  I2 
r  n 


a  ,  n 
n=l 


Z  X"1  (B*E  )  (B*E„)* 


n=l 


n  r  n 


(Eq.  20) 
(Eq.  21) 


3.8  Conventional  beamforming  plus  reference  hydrophone  noise  suppression 
(CB  +  RH) 


The  mean  output  power  in  this  mode  is  (Sect.  B.4  of  App.  B): 


P  =  B*R  B 


1-|B*Q|2  (B*R  B  •  r)"1 


(Eq.  22) 
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With  eigenvector  and  eigenvalues,  this  was  implemented  using: 


B*R  B 


N 

2 

n=l 


B*E 


(Eg.  23) 


N 

B*  Q  =  I  (B*En)  (Q*En)*  .  (Eq.  24) 

n=l 


3.9  "Maximum  Likelihood"  plus  reference-hydrophone  noise  suppression 

(ML  +  RH) 


With  reference  to  Sect.  B.4  of  App.  B: 


_1  _1 

1  -1 

.  -i  -i  -i 

P  =  (B*R  B) 

1+|B*R  Q 

2  {B*R  B(r  -  Q*R  Q)} 

(Eq.  25) 


With  eigen-quantities 


_i  N  .i  |  | 

B*R  B  =  2  k  |B*En|2  ,  (Eq.  26) 


N 

B*R  Q  =  I  An  (B*En)  (Q*En)*  ,  (Eq.  27) 


_i  N  -i  |  I 

Q*R  Q  =  2  X  Q*E  2  -  (Eq.  28) 

n=l  1 


4  RESULTS 

This  chapter  presents  and  discusses  sample  data  from  all  runs  except  Run  2, 
which  are  similar  to  those  from  Runs  1,  3  and  4.  Some  data  are  presented 
from  a  previous  sea-trial  to  show  the  interference  effects  in  the  signal 
paths  from  a  remotely  deployed  transducer  to  an  array. 


4. 1  Signal  path  from  towship  to  array  in  shallow-water 

The  geometry  of  the  trials  is  shown  on  the  left  of  Table  2.  Typical 
geometry  with  an  operational  towed  array  is,  however,  more  likely  to  be  as 
shown  on  the  right  of  the  table. 
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TABLE  2 


GEOMETRY  OF  TRIALS  AND  TYPICAL  GEOMETRY 


Trials 

Typical 

Source  depth 

3  m 

2  to  5  m 

Tow-cable  scope 

200  to  700  m 

500  to  2000  m 

Array  length 

40  m 

50  to  200  m 

Array  depth 

50  to  100  m 

— 

Water  depth 

100  to  500  m 

Together  with  the  type  of  sea-bottom,  the  sea-state,  and  the  sound-speed 
profile,  the  geometry  defines  the  towship's  noise  field  at  an  array  towed 
in  shallow  water.  It  will  thus  be  appreciated  that  the  noise  field  at  the 
towed  array  can  be  very  complex  and  difficult  to  predict.  Thus  care  must 
be  taken  in  using  any  of  the  various  available  techniques. 

Ray- tracing  and  normal -mode  propagation  are  derived  theoretically  in 
App.  D,  indicating  that  ray-tracing  cannot  describe  accurate  low-frequency 
propagation  in  shallow  water.  For  ray-tracing  there  is  a  general  rule  of 
thumb  that  the  water  depth  should  be  greater  than  10A.  Table  3  indicates 
some  safe  limits  for  the  use  of  ray  tracing. 

TABLE  3 

MINIMUM  WATER  DEPTH  (10A)  FOR  SAFE  USE  OF  RAY  TRACING 


Frequency 

Minimum  depth 

(Hz) 

(m) 

50 

300 

100 

150 

250 

60 

500 

30 

Thus,  for  example,  ray  tracing  for  frequencies  of  100  Hz  should  not  be  used 
in  water  shallower  than  150  m;  equally,  in  a  water  depth  of  150  m  ray 
tracing  should  not  be  used  for  frequencies  below  100  Hz. 

In  the  case  of  the  normal -mode  solution  there  is  a  difficulty  at  short 
range  because  of  the  incomplete  description  of  the  acoustical  field  in 
terms  of  the  normal-mode  solution,  which  is  valid  for  farfield  only.  It  is 
usual  to  discard  the  short-range  continuous  field.  Thus  there  is  a  minimum 
range  below  which  the  discrete  modes  alone  do  not  sufficiently  describe  the 
field.  This  minimum  range  is  usually  taken  as  3  to  5  water  depths.  Thus 
in  150  m  water  depth,  ray  tracing  would  be  satisfactory  at  frequencies 
above  100  Hz  out  to  ranges  of  between  450  to  750  m.  Beyond  this  range  the 
normal-mode  solution  should  be  used,  particularly  for  shallow  water.  As, 
in  fact,  different  modes  travel  with  different  velocities  and  arrive 


17 


SACLANTCEN  SM-158 


at  different  angles,  a  mismatch  of  the  array  to  the  acoustical  field  might 
become  a  serious  problem  for  signal  processing  at  low  frequencies.  This  is 
evidently  true  for  the  endfire  direction  of  a  horizontal  array  <5>.  A 
short  explanation  of  the  need  to  distinguish  between  rays  and  modes  for 
array  processing  at  low  frequencies  in  shallow  water  can  be  given  readily 
by  use  of  Eq.  4,  this  being  the  steering  vector  determining  the  beam 
direction.  The  factor  2nxn/A  in  that  equation,  which  is  essentially  the 

array  aperture  size  in  wave  length,  can  be  rewritten  as  u>xn/c  ,  where  c 

is  the  phase  velocity  at  the  receiving  hydrophone  and  is  a  constant  value 
at  one  receiver  for  all  the  rays  impinging  on  this  receiver.  However  for 
all  modes  exhibiting  their  individual  phase  velocities,  this  conventional 
array  processing  will  be  matched  to  one  velocity  only,  hence  with 
increasing  dispersion  the  array  will  become  increasingly  mismatched. 

It  is  not  the  intention  of  this  present  memorandum  to  make  an  extensive 
analysis  of  the  towship's  noise  field  at  the  array  (such  an  analysis  should 
probably  use  the  fast  field  program  or  the  range-dependent  parabolic 
equation  techniques  <4>),  but  merely  to  show  that  the  noise  field  is 
complex,  with  a  covariance  matrix  of  rank  >1.  This  is  the  same  as  saying 
that  the  noise  from  one  source  arrives  from  different  directions  at  the 
same  receiver  when  added  coherently.  The  resolution  of  various  paths  from 
the  towship  with  the  array  was  shown  in  <5>,  using  different  signal- 
processing  techniques.  With  those  methods  and  a  fairly  good  signal-to-noise 
ratio,  only  four  phones  of  a  towed  array  were  used  to  resolve  the  towship's 
noise  directions. 

An  additional  indication  of  the  complex  nature  of  acoustic  signals  in 
shallow  water  was  obtained  as  a  side  result  of  a  previous  sea  trial. 
Figure  7  gives  an  experimental  example  of  mode  interference  over  a  signal 
path  of  3  to  4  km.  The  range  from  the  transmitting  site  to  the  receiving 


~TB  v ai  w  i.  ■  winw  ^  n m  tt.nJl  ■  n  H  E  (j  ■■  ffr*rj  Jl  ■  1 1  ri  ■  M :  F t  ~  - T. i 


xv-  ^  — — 7T   


-  3000  m  - 

10  km  7  km 

DISTANCE  FROM 
TRANSMITTER 


FIG.  7  EXAMPLE  OF  MODE  INTERFERENCE  AT  334  Hz 
RECEIVED  WITH  A  TOWED  ARRAY 
(taken  from  earlier  experiments) 
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ship  varied  from  7  to  about  14  km  during  the  measurements.  For  this 
distance  the  modes  can  be  assumed  to  be  well  separated.  A  strong, 

periodically  repeated  change  (~10  dB)  of  the  received  signal  level  was 
observed  while  the  ship  was  approaching  the  MANNING'S  position,  where  a 
transducer  submerged  at  30  m  water  depth  transmitted  a  334  Hz  CW  signal. 
The  array  was  towed  at  45  m  depth  in  a  shallow-water  area  of  about  80  m 
water  depth  with  a  slightly  upward-sloping  bottom.  The  interference 
wavelength  is  calculated  using  the  approximation  for  a  'soft'  bottom  of 

m  ~  8  H2f/c1(n2-m2).  From  a  comparison  of  the  experimental  and 

theoretical  results  for  the  first  three  modes,  it  is  believed  that  the 
first  and  second  modes  have  the  strongest  influence  on  this  change  in 
level.  Table  4  compares  the  calculated  and  measured  interference 

wavelength. 


TABLE  4 

INTERFERENCE  WAVELENGTH  X  FOR  TWO  DIFFERENT  WATERDEPTHS 
 n,m 


THEORETICALLY 


Mode 

1+2 

2+3 

1+3 

H  =  70  m 

2909  m 

1745  m 

1091  m 

H  =  80  m 

3800  m 

2280  m 

1425  m 

OBSERVED 

H  =  70  m 

H  =  80  m 

~  3000  m 

~  3700  m 

No  further  investigation  was  made  during  the  sea  trial.  The  observed 
values  shown  in  Fig.  7  indicate  the  range  from  own  ship  to  the  transmitting 
site. 

In  the  present  sea-trial  the  water  depth  was  generally  greater  than  150  m, 
the  frequency  over  100  Hz,  and  the  cable  scope  of  the  towed  array  less  than 
750  m.  Therefore  straight-line  ray  tracing  rather  than  normal-mode 
arrivals  was  used  to  estimate  the  angle  of  arrival  at  the  array  of  the 
most  likely  paths  <7>: 

Direct  (D) 

Bottom-bounce  (B) 

Bottom/surface  (BS) 

Bottom/surface/bottom  (BSB) 

Bottom/surface/bottom/surface  (BSBS) 

Figure  8  shows  the  Maximum  Likelihood  Beamformer  (MLB)  output  and  the 
corresponding  water  depths  and  arrival  angles  for  (a)  160  Hz  in  Run  1, 
(b)  335.2  Hz  in  Run  3,  and  (c)  280  Hz  in  Run  4.  The  angles  given  on  the 
sea-depth  profile  correspond  to  the  paths  quoted  above. 
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MAXIMUM-LIKELIHOOD 

BEAMFORMING 

FREQ.  1 60Hz 

TOW  SHIP 

NOISE  TARGET 


0°  30°  I  60°  90°  120°  Il50°  180° 


CABLE  SCOPE  700m 
ARRAY  DEPTH  50m 


TIME  DEPTH  (m) 

o  200  400  600 


FIG .  8a  MAXIMUM  LIKELIHOOD  BEAMFORMING  OUTPUT  SHOWING  CORRESPONDING 
WATER  DEPTHS  AND  ARRIVAL  ANGLES  AT  THE  BOTTOM 
160  Hz  (Run  1) 
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MAXIMUM-LIKELIHOOD 

BEAMFORMING 


FREQ.  335.2Hz 


TOW  SHIP 
NOISE  PATHS 

D  b|  bsb 

0°|  30°  60°f  90° 
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ARRAY  DEPTH  50m 

TIME  DEPTH  (m) 


FIG.  8b  MAXIMUM  LIKELIHOOD  BEAMFORMING  OUTPUT  SHOWING  CORRESPONDING 
WATER  DEPTHS  AND  ARRIVAL  ANGLES  AT  THE  BOTTOM 
335.2  Hz  (Run  3) 
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MAXIMUM-LIKELIHOOD 

BEAMFORMING 

FREQ.  280Hz 


TOW  SHIP 
NOISE  PATHS 

D  B  BSB 


CABLE  SCOPE  500m 
ARRAY  DEPTH  50m 

TIME  DEPTH  Cm) 


FIG .  8c  MAXIMUM  LIKELIHOOD  BEAMFORMING  OUTPUT  SHOWING  CORRESPONDING 
WATER  DEPTHS  AND  ARRIVAL  ANGLES  AT  THE  BOTTOM 
280  Hz  (Run  4) 
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It  is  seen  from  these  figures  that  the  bottom-bounce  path  (B)  is  usually 
the  dominant  one,  although  the  direct  path  (D)  and  the  bottom/surface/ 
bottom  path  (BSB)  are  also  observed,  particularly  at  the  higher  frequency 
when  the  angular  resolution  is  better. 

Figure  8b  includes  data  obtained  whilst  the  towship  was  making  a  90°  turn 
and  it  is  seen  how  the  array  at  500  m  scope  lagged  behind  the  MPG  during 
the  turn  so  that  the  direct  path  made  a  bigger  angle  with  the  array.  These 
direct-path  data  indicate  how,  at  this  frequency,  when  the  array  is  looking 
end-fire  it  can  indicate  the  vertical  angle  of  arrival  of  a  signal.  For 
example,  the  direct  path  under  the  steady-state  conditions  of  Run  3  arrives 
at  about  5°  (arcsin  40/450),  as  shown  in  Fig.  8b. 


4. 2  Towship-noise  suppression 

This  discussion  starts  with  some  mathematical  observations  about  the  three 
methods  implemented: 

a.  Modified  eigenval ues 

b.  Reference  beam 

c.  Reference  hydrophone 

The  question  is  to  determine  how  the  methods  modify  the  covariance  matrix  R 
given  by: 


R  = 


N 

I 

n=l 


A  E  E* 


n  n  n 


(Eq.  29) 


The  method  of  modified  eigenvalues,  as  described  in  Eq.  12,  results  in  a 
modified  matrix  R'  given  by 


R' 


N 

=  1 
n=l 


VE  E* 


n  n  n 


which  can  be  written  as: 


R1  =  R  - 


=  R  - 


(Eq.  30) 


(Eq.  31) 


The  method  with  reference  beam 


(RW  )(RW  )* 

R1  =  R  -  - - —  =  R  -  H2V2V2*  =  R  -  Rt  9  (Eq.  32) 

W*  R  W  1  '6 

r  r 
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The  method  with  reference  hydrophone 


R'  =  R  -  P*  =  R  -  P3V3V3*  =  R  -  Rj  o  •  (Eq.  33) 

y 


Rj  R-j-  ^  and  R-j-  ^  respectively,  are  what  the  three  methods  estimate  as 

the  contribution  to  the  covariance  matrix  R  from  the  towship.  According  to 
definition  (App.  D),  R^  ^  will  have  rank  M  if  M  eigevalues  are  modified, 

while  R-p  2  ar|d  Ry  3  are  seen  always  to  have  rank  1  (only  one  non-zero 

eigenvalue).  But  it  is  shown  in  App.  C  that  the  rank  is  1  only  for  a 
completely  coherent  reception.  Hence  the  two  last  methods  may  work  only  if 
the  real  contribution  R-j-  q  has  rank  1.  The  experimental  results,  however, 

indicate  strongly  that  this  is  not  so. 

The  three  methods  are  discussed  individually  below.  Figure  9  gives  the 
line-printer  outputs  supporting  the  discussion,  although  no  data  are  shown 
for  the  Reference  Hydrophone  Technique  (Sect.  4.2.3). 


4.2.1  The  method  with  modified  eigenvalues 
If  M  eigenvalues  are  modified,  then: 


N 

1 

n=l 


(A  -A'  )  E  E* 
n  n  n  n 


n  £  { n^ ,  n2>  • • •  > 


(Eq.  34) 


Hence  this  method  can  easily  create  a  modification  of  any  rank.  But  that 
is  all;  almost  as  expected,  the  method  is  more  or  less  useless.  It  was 
implemented  simply  because  it  was  so  simple  to  implement.  The  obvious 
reason  that  it  does  not  work  well  is  that  there  is  not  a  one-to-one 
correspondence  between  eigenvalues  vectors  and  sources.  The  contribution 
from  one  source  is  usually  distributed  on  many  eigenvalues,  and  one 
eigenvalue  is  a  combination  of  contributions  from  several  sources.  The 
only  simple  case  is  with  orthogonal  sources. 

Still,  the  method  gave  insight.  It  was  quite  interesting  to  see  the 
"correlation"  between  the  beam-steering  vectors  and  the  eigenvectors: 


cu 


Bi  Ej 


0  “  C-j  j  =  1 


(Eq.  35) 


when  both  B  and  E  are  normalized.  During  the  sea-trial, 


C.  . 


was  printed 


out  for  the  two  eigenvectors  belonging  to  the  two  largest  eigenvalues  (part 
of  the  tabulated  data  given  to  create  Fig.  6).  When  signals  were 
transmitted  from  the  source  below  the  workboat  it  was  often  possible  to 
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CONVENTIONAL  BEAMFORMING 
EFFECT  OF  NOISE  SUPPRESSION  TECHNIQUES 


NO  SUPPRESSION 


SUPPRESSION  BY 
MODIFIED 
EIGEN  VALUES 


SUPPRESSION  BY 
REFERENCE  BEAM 


TOW  SHIP  NOISE 
SUPPRESSION 


TOW  SHIP  NOISE 
SUPPRESSION 


a) 


TOW  SHIP  NOISE 
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320Hz 


a 


TOW  SHIP 

1 80°  0°  90° 


TOW  SHIP  NOISE 
SUPPRESSION 

180*01  90* 


180° 


FIG.  9a  COMPARISON  OF  NOISE  SUPPRESSION  TECHNIQUES 

Run  1  (160  and  320  Hz)  processed  with  conventional  beamforming 
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MAXIMUM-LIKELIHOOD  BEAMFORMING 
EFFECT  OF  NOISE  SUPPRESSION  TECHNIQUES 


NO  SUPPRESSION 


SUPPRESSION  BY 
MODIFIED 
EIGEN  VALUES 


SUPPRESSION  BY 
REFERENCE  BEAM 


FREQ. 

160Hz 


FREQ. 

320Hz 


FIG .  9b  COMPARISON  OF  NOISE  SUPPRESSION  TECHNIQUES 

Run  1  (160  and  320  Hz)  processed  with  maximum  likelihood  beamforming 
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CONVENTIONAL  BEAMFORMING 
EFFECT  OF  NOISE  SUPPRESSION  TECHNIQUES 


NO  SUPPRESSION 


SUPPRESSION  BY 
MODIFIED 
EIGEN  VALUES 


SUPPRESSION  BY 
REFERENCE  BEAM 


TOW  SHIP  NOISE 


TOW  SHIP  NOISE 
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c) 
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FIG .  9c  COMPARISON  OF  NOISE  SUPPRESSION  TECHNIQUES 

Run  2  (280  and  334.4  Hz)  processed  with  conventional  beamforming 
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MAXIMUM-LIKELIHOOD  BEAMFORMING 
EFFECT  OF  NOISE  SUPPRESSION  TECHNIQUES 


NO  SUPPRESSION 


SUPPRESSION  BY 
MODIFIED 
EIGEN  VALUES 


SUPPRESSION  BY 
REFERENCE  BEAM 


TO W  SHIP  NOISE 


TOW  SHIP  NOISE 
SUPPRESSION 

30”/  60°  90“  120°  150°  180“ 


FIG.  9d  COMPARISON  OF  NOISE  SUPPRESSION  TECHNIQUES 

Run  2  (280  and  334.4  Hz)  processed  with  maximum  likelihood  beamforming 
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observe  how  the  two  eigenvalues  Ex  and  E2  were  correlated  with  beam 

directions  corresponding  both  to  the  towship  noise  and  to  the  workboat 
source.  This  means  that  we  could  easily  suppress  the  towship  noise  by 

modifying  the  first  three  to  four  eigenvalues,  but  at  the  same  time  the 

other  signals  would  also  be  suppressed  in  an  uncontrolled,  random,  and 

untolerable  way. 

At  this  point  it  may  be  appropriate  to  speculate  why  the  towship  noise  seen 
by  the  array  towed  500  m  behind  the  ship  has  a  covariance  matrix  of  rank  3 
to  4.  There  are  three  possibilities: 

a)  That  the  geometry  of  the  array  relative  to  the  ship  is  changing 

during  the  averaging  time  of  the  covariance  matrix.  This  time  was  40  to 
60  s  in  the  sea-trials,  which  is  obviously  long  compared  with  the  possible 
snaking  of  the  array.  The  solution  would  be  to  use  a  much  shorter 
averaging  time  and,  if  possible,  to  overlap  the  FFT-windows  in  the 
filtering  process  before  the  averaging,  so  as  to  use  the  data  efficiently. 

b)  That  the  different  paths  of  the  multipath  propagation  from  the 
towship  to  the  array  have  different  travel  times.  As  the  continuous 
background  spectrum  of  the  towship  has  a  short  correlation  time,  the 
contribution  at  a  given  time  from  the  different  paths  may  have  a  reduced 
cross-correlation,  which  would  mean  a  non-coherent  signal  from  the  towship 
and  thus  a  covariance  matrix  of  rank  higher  than  one. 

c)  That  the  influence  of  the  medium  and  its  boundaries  is  reducing 
the  coherence,  even  for  an  originally  stable  line  (CW). 


4.2.2  The  method  with  reference  beam 

This  method  estimates  the  covariance  matrix  of  towship  noise  to  be 


(RWr)(RWr)* 

By  2  ~  W*~RW  *  ^2^2 ^2*  »  (Eq.  36) 

*  r  r 

where 


W 


r 


*  _i  _i  _i 

(BrR  Br)  R  Br 


for  CBF 
for  MLBF  . 


In  both  cases  the  model  assumes  that  the  towship  noise  arrives  as  a 
coherent  plane  wave  at  the  array.  This  might  work  well  in  deep  water;  in 
shallow  water  it  definitely  did  not.  The  towship  noise  was  certainly 
suppressed  in  a  bearing  interval  around  the  reference-beam  direction  (given 
by  B^) ,  but  the  rest  of  the  towship  noise  was  more  or  less  unchanged. 

Secondly,  if  there  was  any  signal  (target)  in  the  reference  direction,  it 
would  also  be  cancelled.  The  plots  indicate  multipath  arrival  at  different 
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vertical  arrival  angles.  Hence,  even  if  the  towship  noise  is  coherent,  it 
would  not  fit  a  plane-wave  model.  Making  it  even  more  difficult,  there  are 
strong  indications  (from  study  of  eigenvalues  and  vectors)  that  the 
towship-noise  signal  is  not  coherent  over  the  array. 

It  might  be  possible  to  work  out  an  ad  hoc  algorithm  by  choosing  more  than 
one  reference  beam  and  cross  correlating  between  the  reference  beams,  but 
it  is  not  obvious  that  the  outcome  would  be  worth  the  effort. 


4.2.3  The  method  with  reference  hydrophone 

This  method  estimates  the  covariance  matrix  of  the  towship  noise  to  be 


where 


RT,3  =  M.V.V 


Q  =  ^  I  X.y.= 
M  .=1  r*! 


M 


and 


(Eg.  37) 


(Eq.  38) 


As  seen,  the  method  assumes  the  covariance  to  be  of  rank  1,  which  means  a 
coherent  signal.  But  it  is  an  improvement  (theoretically,  at  least)  from 
the  reference-beam  method  in  that  it  does  not  assume  anything  about  the 
shape  of  the  wavefront.  Furthermore,  even  a  signal  with  the  same 
directivity  pattern  as  the  towship  noise  should  be  only  slightly  reduced 
if: 


ct2  /a2  \ 

ref/  array J Towship 


» 


o2  /  a2  \ 
ref/  array  / 


Target  , 


where 


°2re^  =  Towship  or  Target  power  on  the  reference  hydrophone 

°2 array  =  Towsh'iP  or  Target  power  on  the  array  hydrophone. 

This  is  exactly  the  reason  why  the  reference  hydrophone  has  to  be  closer  to 
the  towship  than  to  the  array.  From  a  theoretical  point  of  view  this 
should  be  the  best  of  the  three  implemented  methods.  Although  this  sounds 
promising,  the  results  so  far  show  that  the  method  in  its  present  shape 
does  not  work  at  all.  There  are  two  possible  reasons  for  this: 

a)  The  first  is  that  the  reference  hydrophone  was  rather  simply 
designed  and  thus  created  flow  noise.  This  would  especially  influence  its 
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power  estimate  r^  .  This  estimate  was  therefore  modified  by  a  factor 
smaller  than  1,  to  see  what  would  happen: 


r  '  =  q  r  ,  where  OS  q  SI. 
y  n  y’ 

This  was  done  at  different  frequencies,  including  the  frequency  of  334  Hz 
transmitted  by  the  workboat. 

The  result  was  disastrous.  The  workboat  signal  disappeared  almost 
completely,  while  the  towship  noise  was  still  there,  reduced  only  for  some 
bearings.  Looking  a  little  closer  at  the  mathematics,  let  By  be  a 

direction  where  we  can  identify  towship  noise,  and  B,.  be  the  direction  of 

the  workboat  source.  Then,  in  obvious  notation: 


Py  =  B*  R  By 


-11M„  _  n*  „  n  _  l“T 


B*  Q|2 

PT  =  BT  lR  "  V-  ]BT  =  BT  R  B  "  r 


P* 


1  |BfQ|2 


r  B*  RBT 
y  T  T 


Similarly: 


Pi 


_S  _y  1  |BS  Q|2 
PS  ry  BS  R  BS 


(Eq.  39) 
(Eq.  40) 


(Eq.  41) 


(Eq.  42) 


To  obtain  a  good  suppression  of 
suppression  of  the  workboat  signal,  we 

PT 

<< 

P‘ 

_s 

PS 

which 

is  the 

same 

as: 

B*  Q|2 

j 

Bs  ^l2 

B*  R  By 

B*  R  Bs 

the  towship  noise  relative  to  the 
must  have: 


(Eq.  43) 


(Eq.  44) 


This  is  in  many  (maybe  most)  situations,  not  the  case. 

The  second  possibility  is  that  the  change  of  the  geometry  relevant  to 
towship  noise  (towship  source,  reference  hydrophone,  array)  during  the 
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averaging  time  may  have  been  much  more  damaging  than  the  changes  in 
geometry  relevant  to  the  workboat  signal.  There  would  always  have  been 
some  heaving  movement  of  the  towship  that  could  have  changed  the  position 
of  the  noise  source  and  the  nature  of  the  multipath  pattern  of  its  noise 
propagation.  The  implemented  averaging  time  was  certainly  too  long  to 
allow  an  adjustment  to  the  dynamics. 

Even  if  the  geometry  was  completely  frozen,  the  present  implementation 
might  still  not  work  well.  The  reason  is,  as  mentioned  earlier,  that  the 
rank  of  the  towship  noise  at  the  array  was  higher  than  1,  at  least  of  the 
order  3  to  4  at  334  Hz.  Hence  it  cannot  be  suppressed  by  a  rank-1  matrix. 

There  are  two  ways  to  increase  the  rank  of  the  reference  system.  One  is  to 
increase  the  number  of  reference  hydrophones.  This  might  help  to  some 
extent,  but  it  will  not  help  if  the  reason  for  reduced  coherence  is  the 
different  travel  times  along  the  different  paths.  Another  possibility  is 
to  work  with  several  delayed  versions  of  the  reference  hydrophone  output, 
to  create  a  vector: 


(  y(t-tx) 

y(t)  =  <  y(t-t2)  .  (Eq.  45) 

(  y(t-tL) 


By  going  through  the  mathematics,  this  will  result  in: 


R'  =  R 


xx 


R  R 

*y  yy 


_i 


R  =  R-R 
yx 


T,y 


where 


1  M 

R  =  jjj  I  X .  Y* 
xy  M  i=1  ^  i 


etc. 


(Eq.  46) 


(Eq.  47) 


The  dimension  L  of  Y  should  be  such  that  Ryy  has  full  rank.  Otherwise  it 

cannot  be  inverted,  and  in  fact  has  more  elements  than  can  be  used.  Ideally 
the  rank  of  Ryy  should  be  the  same  as  the  rank  of  the  towship-noise 

contribution  to  R. 

The  problems  of  flow  noise  could  be  overcome  by  creating  a  new 
reference-sensor  system  with  reduced  flow  noise  and  low,  or  almost  no, 
sensitivity  to  acceleration  (vibration). 


4. 3  Target  indication  and  bearing  estimation 

As  with  the  techniques  for  the  suppression  of  towship  noise,  only 
qualitative  comments  can  be  given  on  the  performance  (Fig.  10)  of  the  two 
forms  of  processing  (CBF,  MLB)  for  providing  an  indication  of  the  presence 
of  a  target  of  opportunity  and  an  estimate  of  the  target  bearing  for 
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frequencies  of  112  Hz  and  334.4  Hz.  Figure  10b  also  includes  the  presence 
of  a  334.4  CW  source.  Points  to  note  from  Figs.  10a  and  b  are: 

a.  The  CBF  processing  does  not  include  any  shading,  therefore  the 
sidelobes  are  high.  However,  CBF  is  not  a  serious  contender  for 
providing  target  data. 

b.  The  MLB  beamforming  gives  a  sharp  and  precise  display,  although 
at  times  low-level  'targets'  are  suppressed,  due  to  the  maximum  level 
used  on  the  normalizing  factor  in  the  display. 


CONCLUSIONS 

a.  A  detailed  analysis  should  be  made  of  the  signal  paths  between 
the  towship  and  the  array.  An  eigenvalue  analysis  indicated  that  the 
multipath  propagation  is  reducing  the  coherence  (Sect.  4.2).  This  may 
require  the  application  of  the  FFP  technique  and  normal -mode  propagation 
(Sect.  4.1). 

b.  Because  of  the  lack  of  a  one-to-one  correspondence  between 
eigenvalue/vector  and  source,  except  when  the  sources  are  orthogonal,  the 
modified  eigenvalue  technique  for  suppressing  towship  noise  does  not  appear 
to  be  feasible  (Ch.  3).  However,  the  eigenvalues/vectors  give  useful 
information  on  certain  characteristics  of  signals  and  will  therefore  be 
used  in  the  future  signal-processing  software. 

c.  The  use  of  a  single  reference  beam  for  suppressing  towship  noise 
is  not  satisfactory;  it  assumes  that  the  towship  noise  arrives  at  the  array 
as  a  coherent  plane  wave  and  that  there  is  no  target  signal  in  that 
direction  (Ch.  3).  It  would  be  possible  to  select  more  than  one  beam  and 
produce  an  ad  hoc  algorithm  but  its  usefulness  under  the  full  range  of 
operating  conditions  would  be  questionable.  No  further  work  is  proposed  at 
this  stage. 

d.  The  reference-hydrophone  technique  did  not  work  because  of  the 
oversimplified  nature  of  the  present  technique  and  the  high  level  of 
flow-noise  in  the  reference  hydrophone  (Ch.  3).  A  technique  to  provide  a 
more  accurate  representation  of  the  towship  noise  is  being  developed  for 
future  sea  trials,  and  the  hydrophone  design  is  being  improved  to  keep  the 
flow-noise  level  to  a  minimum  using  a  separate  streamer  for  the  reference 
hydrophone. 
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TS  TOW  SHIP  FREQUENCY  112Hz 

T  TARGET 

CONVENTIONAL  MAXIMUM-LIKELIHOOD 

BEAMFORMING  BEAMFORMING 


FIG .  10a  COMPARISON  OF  TARGET-BEARING  INDICATION  BY  CONVENTIONAL 
BEAMFORMING  AND  MAXIMUM  LIKELIHOOD  BEAMFORMING 
Low  frequency ,  112  Hz,  Target  only 
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TS  tow  SHIP 
T>  TARGET 

S  source  FREQUENCY  334.4Hz 

SL  SOURCE  LEVEL  CdBl 


CONVENTIONAL  MAXIMUM-LIKELIHOOD 

BEAMFORMING  BEAMFORMING 


FIG .  10b  COMPARISON  OF  TARGET -BEARING  INDICATION  BY  CONVENTIONAL 
BEAMFORMING  AND  MAXIMUM  LIKELIHOOD  BEAMFORMING 
High  frequency,  334.4  Hz,  Target  plus  CW  source 
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APPENDIX  A 

COVARIANCE  MATRIX  ESTIMATE 


This  appendix  gives  seven  definitions  with  their  associated  properties  and 
"proofs"  relevant  to  our  use  of  the  covariance  matrix  estimate,  R.  The 
main  emphasis  is  on  eigenvalues  and  eigenvectors. 

The  square  matrix  R  of  dimension  N  is  defined  by: 


M 

2  X.  X? 
,  1  1 


(Eq.  A. 1) 


where  X  will  be  called  a  data  vector,  and  is  of  dimension  N,  and  the  star 
indicates  a  complex  conjugate  transpose. 


A.  1  DEFINITION 

A  square  matrix  P  is  called  hermitian  if  it  is  equal  to  its  own  complex 
conjugate  transpose,  P*. 

Property  1:  R  is  a  hermitian  matrix.  (PI) 


Proof: 


/  M 

I  X.  X5!? 
\i=l  1  1 


M 

2 

i=l 


X. 


(Eq.  A. 2) 


Property  2  Z  =  Y*RY  is  a  real  scalar  for  all  vectors  Y:  (P2) 

Proof: 

1*  =  (Y*RY)*  =  Y*R*  Y 
By  use  of  (PI): 

Z*  =  Y*!R* Y  =  Y*RY  =  Z  q.e.d. 


A. 2  DEFINITION 

A  square  hermitian  matrix  P  is  positive  definite  if,  for  all  vectors 
Y  except  the  null  vector,  Y*PY  >  0,  and  non-negative  definite  if  Y*PY^  0. 
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Property  3:  R  is  a  non-negative  definite  matrix. 


(P3) 


Proof: 


(Eq.  A. 3) 


A. 3  DEFINITION 

If  E  is  a  vector,  but  not  a  null  vector,  such  that  PE  =  AE,  then  E  is 
called  an  eigenvector,  and  A  the  corresponding  eigenvalue.  The  eigen¬ 
vectors  are  always  assumed  to  be  normalized:  E*E  =  1. 


Property  4:  the  eigenvalues  of  R  are  real  and  non-negative. 


(P4) 


Proof: 

AE  =  RE 

Premultiplying  with  E*: 

A  =  A  E*  E  =  E*  R  E 

According  to  (P2),  A  is  real,  and  according  to  (P3),  A  is  non-negative. 

A. 4  DEFINITION 

Two  vectors  X  and  Y  ,  not  being  nul  1- vectors ,  are  orthogonal  if  X*Y  =  0. 
They  are  called  orthonormal  if  both  X  and  Y  are  normalized  by  having  unit 
length. 

Property  5:  Eigenvectors  of  R  with  different  eigenvalues  are 

orthogonal . 

Proof: 


By  exchanging  indices: 
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Then  take  the  complex  conjugate  transpose  of  both  sides  of  the  last 
equation  and  use  (PI)  and  (P4): 


A  E*E  =  E*  R  E 
n  n  m  n  m 


Combining  the  first  and  last  equations  of  the  proof: 


Xm  - 


A  )  E* 
n  '  n 


Em  = 


As  the  eigenvalues  are  assumed  different  A  /  A  : 

M  n  m 


E  * 
n 


E  =  0 
m 


q.  e.  d. 


Property  6:  There  are  N  eigenvalues  of  R  ,  which  may  not  be  all 

distinct  (different  from  each  other),  and  they  are  a  solution  to 

det(R-AJ)  =  0 


Proof: 

RE  =  AE 
(R-AI)  E  =  0 

This  can  be  true  for  E  not  being  a  null  vector  if,  and  only  if,  the 
determinant  of  (R-AI)  is  zero,det(R-AI)  =  0. 

This  is  a  polynome  in  A  of  degree  N  ,  and  has  N  solutions,  which  may 
not  be  distinct.  These  N  solutions  are  hence  eigenvalues  of  R. 

We  know  already  that  if  eigenvalues  are  different,  their  corresponding 
eigenvectors  are  orthonormal.  If  there  are  k  eigenvalues  with  the  same 
value,  we  will  state  without  proof  that  it  is  always  possible  to  find 
exactly  k  orthonormal  vectors  that  are  eigenvectors  of  R  with  this  same 
value  on  their  eigenvalues.  Hence  we  state  without  proof  the  following 
property: 


Property  7:  R  has  always  N  orthonormal  eigenvectors.  (P7) 

Property  8:  The  eigenvectors  with  positive  eigenvalues  are  linear 

combinations  of  the  data  vectors. 

Proof: 

,  M  ,  M 

AE  =  RE  =  rj  I  X.X.E  =  n  *  (X.E)X.  q.e.d. 

n  i=l  11  i=l  1  1 
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Property  9:  Eigenvectors,  which  are 

vectors,  have  positive  eigenvalues. 


linear  combinations  of  the  data 

(P9) 


Proof: 


,  M  M 

K  =  A.E*E  =  E*RE  =  -k  2  E*X.X.E  =  ^  2 

M  i=i  '  1  M  f=i 


X*E 


>  0 


q.  e.  d. 


A. 5  DEFINITION 

The  k  vectors  Y  are  said  to  be  linearly  independent  if 
K 

2  c.Y.  =  0  only  if  c.=  0  for  i  =1,  ...,  k 

i=l  1  1  1 

Property  10:  If  the  number  of  independent  data  vectors  is  M,  and  M  <  N, 

N  -  M  eigenvectors  are  orthogonal  to  the  data  vectors,  and  their  corres¬ 
ponding  eigenvalues  are  zero.  (P10) 

"Proof" : 

From  M  independent  data  vector,  we  state  without  proof  that  M  eigen¬ 
vectors  be  constructed: 

M 

E  =  I  c.X,  m  =  1,  2,  ... ,  M 

m  i=1  i  i 


or,  in  a  more  compact  form: 

[El*  E2’  •••»  Em]  =  [  Xl»  X2’  **•*  Xm]  C 


As  the^  M  eigenvectors  are  also  independent,  C  must  be  non-singular, 
and  C  exists: 

[Xl’  X2 . X  m]  [  El-  E2 . En]  C  1 

Premultiply  with  an  eigenvector  En  n  =  M+l,  ...,  N 

[ESX1>  EnX2 . ES  Xn]  =  [En  Er  ES  E2 . ES  Em]  C_1 

=  [o,  0,  ...,  o]  c'1  =  0 


by  using  (P7).  Hence 

E*Xi  =0  for  i  =  1,  2,  ... ,  M 


n  =  M+l,  M+2,  . . . ,  N  q.e.d. 
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The  corresponding  eigenvalues  are: 


=  A  E*E 


n  n  n 


E*R  E 
n  n 


q.  e.  d. 


A. 6  DEFINITION 

A  square  matrix  P  is  singular  if  its  determinant  is  zero,  and  its  inverse 
P  1  does  not  exist. 


Property  11:  If  the  number  of  independent  data  vectors  M  is  less 

than  N,  R  is  singular.  (Pll) 


Proof:  From  (P10)  it  follows  that  if  M  <  N,  N-M  eigenvalues  are 

zero.  But  the  eigenvalues  are  solutions  to  det(R-Al)  =  0. 

Hence  for  A  =  0  ,  det(R)  =  0  and  R  is  singular. 

Property  12:  If  R  has  N  independent  data  vectors,  it  is  non-singular 
and  its  inverse  R  ^  does  exist.  (P12) 

Proof:  The  proof  follows  directly  from  (P8)  and  (P9)  for  M  =  N  and 

these  data  vectors  are  independent. 


A.  7  DEFINITION 

A  square  non-singular  matrix  P  is  Unitarian  if  its  inverse  P  ^  is  equal 
to  its  complex  conjugate  transpose  P*. 

Property  13:  The  square  matrix  S  =  [E^,  E^,  ...,  En]  constructed  from 

the  N  eigenvectors  of  R,  are  Unitarian.  (P13) 

Proof:  Because  the  eigenvectors  are  independent,  S  is  non¬ 

singular,  and  its  inverse  S  ^  exists.  Let  us  introduce: 

T  =  (S-1)* 

=  [Tl>  t2,  ... ,  tn] 

where  T.  is  a  vector. 

These  vectors  can  be  expressed  as  linear  functions  of  the  eigenvectors: 

N 
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From  T*S  =  S  =  I  it  follows  that: 


N 

T .*E  =  I  a* .  E?  E  =  a*.  =1  if  i  =  n 

in  ij  j  n  in 

j=l 

=0  if  i  ^  n 


From  thi s  fol lows: 


-1, 


Ti  =  E.  and  T  =  S  =  (S  )*  q.e 


Property  14:  R  can  be  written  as  R  =  SD  S*,  where 

D  =diag(\1,  A2,  . \  )  is  a  diagonal  matrix. 

Proof: 

RE.  =  E.  i  =  1,  2,  N  can  be  combined  to: 

RS  =  SD 


As  S  is  non-singular: 
R  =  S  D  S"1 


and  from  (P13): 

R  =  S  D  S*  q.e.d. 


Property  15:  If  R  is  non-singular,  its  inverse  can  be  written: 

_  i  _i 

R  =  S  D  S*  . 


Proof: 

R_1  =  (S  D  S*)~ 
Using  (P13),  this  is 
r"1  =  S  D"1  S* 


_i  _i 

(S*)  D 


_i 

S 


q.  e.  d. 


(P14) 


(P15) 
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Property  16:  Let  n  be  an  integer,  then 

n  n  * 

Rn  =  S  Dn  S 

Proof: 


For  n  ^  1 

Rn  =  (SDS*)  (SDS*)  ...  (SDS*)  =  SDS_1  SDS_1  ...  SOS*  =  SDnS* 

For  n  S  -1 

R~n  =  (SD‘1S*)(SD“1S*)...(SD"1S*)  =  SD~VXS  D“V\ .  .SD_1  =  SD 
For  n  =  0 

R°  =  I  =  SS_1=  SS*  =  SD°S*  q.e.d. 


Property  17:  Let 

n 

be  an  integer,  then 

N 

n 

R  =  I  X.  E. 
i=l  1  1 

E* 

i 

Proof: 

Dn  =  [diagCXj, 

^2  > 

•  ••»  ■X.j.j) 3  =  diagCAj,  X2,  •••»  X^) 

by  direct  use  of  matrix  multiplication.  Hence: 


Rn  =  [Ei,  E2,  ....  EJ 


[^1)  ^2  y  •••»  E|^j] 


X? 


n  ^ 

xf  £1 

n  * 

X?  Eo 


XN  EN 


N 

I 

i=l 


* 

Ei 


Xi  Ei  Ei 


q.e.d. 


q.e.d. 


nS*  q.e.d. 
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Property 

X* 

Proof: 


18:  For  all  X  and  Y: 

Rn  Y  =  I  A1?  (X*  E^E*  Y) 
i=l 

This  follows  directly  from  (P17). 
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APPENDIX  B 

BEAMFORMING  AND  NOISE  SUPPRESSION 


B.l  CONVENTIONAL  BEAMFORMING  (CB) 

Let  X.  be  a  narrowband  complex  data  vector  of  dimension  N  (N  array 

hydrophones)  from  a  frequency  band  around  f,  where  i  is  a  time  index. 
Then  the  covariance  matrix  estimate  R  is  defined  by: 

i  M  * 

R  =  jjj  2  X.  X.  M  £  N  ,  (Eq.  B.l) 

"  i=l  1 

where  1/M  is  added  for  convenience. 

A  steering  vector  B  is  defined  by: 


cos  p 


where  x^  is  the  position  of  the  n-th  hydrophone,  p  is  the  bearing  angle 

measured  from  ahead  endfire,  and  /"FT  is  introduced  to  normalize  B: 

B*B  =  1. 

A  "conventional"  beam  is  then: 

* 

Zk  i  =  B^  X*  ,  where  i  and  k  are  time  and  bearing  indices. 


The  power  output  is: 


P 


k,i 


X.  X*  B. 
i  l  k 


The  mean  (in  time)  conventional  beam  power  output  is: 

M  *  M  *  * 

Pk  =  H  ^  pk,i  '  B  Bk  "i  xi  Bk  *  Bk  R  Bk  • 

Using  eigenvalues  and  eigenvectors,  this  can  be  written: 


Pk  =  j  VBk  En)(En  Bk>  =  \  \ 


* 

Bk  En 
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B. 2  MAXIMUM  LIKELIHOOD  BEAMFORMING  (MLB) 

Instead  of  the  steering  vector  B  ,  use  a  weight  vector  W  to  form  a  beam: 

* 


The  so-called  maximum-likelihood  beamforming  is  obtained  by  choosing  a 
weight  vector  W  that  is  minimizing 

R  constrained  by  Bk  =  B^  =  1 

Y.  =  W*  R  W.  =  (W*  -  W  )  R  (W.  -  W  )  +  wf  R  W  +  W*  R  W.  -  W*  R  W 
k  k  k  k  o  k  o  k  o  o  k  o  o 

The  constraint  can  now  be  introduced  by  choosing 

Yk  =  (WR  -  c  R-1Bk)*  R(Wk  -  c  R_1Bk)  +  c  +  c*  -  cc*  B*  R"^ 


Assuming  R  to  be  positive  definite,  this  has  a  minimum  for: 

W  =  c  R-1  B. 
k  k 

The  constant  c  is  determined  from  the  constraint 

*  *  _  i 

B.  W.  =  c  B.  R  B  =  1 
k  k  k  k 

*  _  i  _  i  _  i 

W .  =  (B.  R  B.  )  R  B. 
k  v  k  k'  k 

The  ML  beamformer  output  will  be: 

*  _i  _i  *  _i 

zm  =  <Bk  R  V  <Bk  R  V 

and  the  mean  power  output  will  be: 

Pk  =  Wk  R  Wk  =  (Bk  R_1  Bk)_1  (Bk  R_1  R  (Bk  R"lBk) 

*  _  1  ,1 

=  <BkR  V 


Using  eigenvectors  and  eigenvalues  of  R: 


N  i  *  a 

I  \  (B.  E  )  (E  B.  ) 

,  n  v  k  n  v  n  k' 
n=i 
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B. 3  NOISE  SUPPRESSION  WITH  REFERENCE  BEAM 


Let  there  be  a  reference  beam: 


V.  =  W  X.  where  W  =  B  for  CBF 

i  r  i  r  r 

=  (B  R~  B  )"  R"  B  for  MLBF 
r  r  r 


This  reference  beam  is  used  to  correct  the  beams: 


Zk,i  =  Zk,i  -  \  Yi 


Vk  is  now  chosen  to  minimize: 


p,'  =  i  i 


k  Mi=l 


M  , 
Z 


k,i 


1  M 

2  =  -  I  (1 
M  k«i 


V.  V.)  (Z.  .  -  V.  Y.) 
k  i  k,i  k  i 


*  * 


M 


Tj  (Zk,i  Zk,i  '  Zk,i  Yi  Vk  "  Vk  Yi  Zk,i  +  Vk  Vk  Yi  V 


=  r2  '  ’’zy  Vk  -  Vk  rzy  +  Vk  \  ry 


with  obvious  notation.  This  can  be  written: 


p.  =(V  -  ^Vv,-  ^Vr  +  r-  -zy  V 
k  ^ k  ry  /  v k  ry  /  y  2  ry 


This  has  a  minimum  for: 


_  ^  %  « vy 

k  r  W*  RW 

y  r  r 


where  Wk  is  defined  similar  to  W  for  CBF  and  MLBF. 


The  corrected  linear  and  mean  power  outputs  are: 


WkRWr  *  /  *  WkRWr  * 

Zk  i  =  Wk  Xi"  JL~1  Wr  Xi  =  Wk"  jL-r  Wr  1  Xi 
k>  k  W*RW  r  1  l  k  W*RW  r'  1 
r  r  \  r  r 


1  *  /  *  \_i  I  * 

k  =  WkRWk-^rRWrj  |wk 


■  R  W 
k  r 
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for  CBF: 


Pk  =  Bk  R  Bk  '  <Br  R  Br^" 


B.  R  B 
k  r 


Using  eigenvalues  and  eigenvectors: 


^  N  -  *  | 

BkRBk  =  Z/„  BkEn2 

n=l 


B  R  B  =  I  X 
r  r  n 

n=l 


N  * 

B  E 


r  n 


B*  R  B 
k  r 


N  a-  *j k 

1  X  (B.E  )(B  E  ) 
,  n  K  k  r  n 

n=l 


For  MLBF: 


i  *  _i  _i 


_i  _i 


.  1  _  2 


Pk  =  <Bk  R  V  -  <Br  R  V  (Bk  R  V 
And  with  eigenvalues  and  eigenvectors: 


*  _i 

B.  R  B 
k  r 


*  - 1  N  _  i 

B.  R  B.  =  X  \ 
k  k  ,  n 

n=l 


* 

B,  E 
k  n 


*  _  x  N  _  i  |  * 

BR  B  =  I  \  BE 
r  r  n  r  n 

n=l  1 


*  _i 

B.  R  B  = 
k  r 


N  _i  /  *  \  /  *  V 

1  k  B.  E  )  B  E  ) 

n=l  n  V  k  "A  r  7 


B. 4  NOISE  SUPPRESSION  WITH  REFERENCE  HYDROPHONE 

Let  there  be  a  reference  hydrophone  with  output  y.  This 
hydrophone  output  is  used  to  correct  the  array  hydrophone  outputs: 

x'.  =  X.  -  Wy. 


W  is  now  chosen  to  minimize: 
M 


1  1  _ 

P  =  „  1 

M  i=l 


x.  -  Wy. 


2  - 


->  M  * 

s  2  -  wyp  (xi  -  wYi) 


M 


i=l 


reference 
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=  -str 


=  ^tr 


M 

I  (X.  -  Wy.)(X.  -  Wy.) 
1=1  1  1  1  1 


M  •OU  «JU  <JL*  X  X  X 

I  (X.  X.  -  X.  y.  W  -Wy.X.+WW  y.y.) 

ii  l  ^i  i  Jt 

i=l 


=  tr  (R  -  QW*  -  WQ*  +  WW*  r  ) 
with  obvious’  notation.  This  can  be  written: 


tr 


K)(-y 


r  +R-2S! 
y  ryj 


This  has  a  minimum  for: 


W 


=  2_ 


The  corrected  data  vector  and  covariance  matrix  will  be: 

x‘.  =  X.  -  2-  y. 
i  i  r  ■’i 

y 


1  1  1  *  nn* 

R1  =  4  I  X.  X.  =  R  - 
M  i  i  r 

1=1  y 


The  new  inverse  matrix  will  be  using  the  matrix  inversion  lemma: 
R' 


M"1  =  ^  =  R-1  +  ( 


ry  -  Q*R"1Q)'1(R"1  QQ*  r"1) 


The  mean  output  power  using  CBF: 


Pk  BkR  Bk  '  BkRBk  '  ry 


* 

B,Q 
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Using  eigenvalues  and  eigenvectors: 


*  ^1*1 
BkRBk%^Xn|Bk  En|2 


o  * 


Bk«  =  (Bk  En>  5  En> 

n-1 


This  way  of  writing  B*Q,  as  B*IQ  =  B*R°Q  may  seem  artificial,  but  may  be 
useful  from  a  software  point  of  view. 

The  mean  output  power  using  MLBF: 


*  _i  _  i 


[*  _  i 

Bk  R 


Pk  -  (Bk  «■  Bk)  ,-k  ••  -k 

Using  eigenvectors  and  eigenvalues: 

*  _  i  N  _  1 1  * 

BkR  K  E„ 


Bk  +  (ry  -  Q^'q)"1  |B*  R 


-,ql 


_  i 


N 


Q*  R  Q  =  I  \r 
n=l 


_  i 


Q*  E  2 

n  I 


*  _i 


BkR  Q  = 


N  x  *  ^ 

t,  \  (Bk  En>  En> 

n=l 
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APPENDIX  C 

THE  RANK  OF  A  SOURCE'S  CONTRIBUTION  TO  THE  COVARIANCE  MATRIX 


From  Appendix  A,  we  know  that  the  covariance  matrix  R  can  be  written: 

N  * 

R  =  SAS*  =  Z  A.  E.  E. 

i=l  1  1  1 


If  only  the  first  M  eigenvalues  are  non-zero,  R  is  singular,  and  is 
said  to  have  rank  M: 

M  * 

R  =  I  A,  E,  E. 


i=l 


i  i  i 


Let  us  look  at  the  contribution  from  one  source, 
completely  coherent  over  the  total  array,  then: 

i  iOJj.-iJj.) 

E{x.x*}  =  (E{|x.|2}  •  Edx^})'*  e  1  J 


If  its  contribution  is 


* 

V.V. 
1  J 


leading  to: 

Rs  =  E{XX*}  =  VV* 

where  V  is  a  normalized  vector  (V*V  =  1).  Now: 

R  E  =  o2  V(V*E)  I  =  °s  V  !f  E*V  -  1  (E  =  V) 
s  s  '1=0  if  E*V  =  0  (E  X  V) 


This  means  that  there  is  only  one  non-zero  eigenvalue  A.  =  with 

corresponding  eigenvector  E,  =  V  and  the  rank  of  R  is  1. 

l  s 

If  the  contribution  from  the  source  is  completely  incoherent,  then: 


E{x  x*}  i 

J  1  =  0 

if 

i  =  j 

if 

i  *  J 

Then: 

Rs  =  diag  {of  .. 

.  o? 

•  ‘  *aN 

As  det  R<.  -  n  a?  >  0  ,  this  matrix  has  always  full  rank,  or  rank  equal  to  N. 
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Hence  a  real  source  will  give  a  contribution  to  the  matrix  with  a  rank  from 
1  in  the  case  of  complete  coherence,  and  increasing  rank  with  decreasing 
coherence.  We  will  define  the  field  from  the  source  to  have  rank  M  if: 


10  lg 


^  a 

for 

i  S  M 

<  a 

for 

i  >  M 

where  a  is  a  chosen  threshold  (as  an  example,  -10  dB)  and  the  eigenvalues 
are  arranged  in  decreasing  order. 
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APPENDIX  D 

RAY  THEORY  AND  NORMAL-MODE  THEORY 


The  acoustic  propagation  is  usually  expressed  in  terms  of  normal  modes  for 
the  low-frequency  propagation  or  in  terms  of  rays  for  the  high-frequency 
range.  Ray  theory  and  normal-mode  theory  is  briefly  reviewed  here.  A 
numerical  comparison  and  overview  is  in  reference  <6>. 


D.l  RAY  THEORY 

The  limiting  form  of  the  exact  wave  theory  for  A  -»  0,  i.e.  f  »  is 
called  ray  theory  and  the  propagation  is  found  to  be  frequency- independent 
for  this  limit.  To  achieve  the  limit  the  wave  equation  for  the  sound 
pressure 


V2p  -  c(z)"2  |E  =  0  (Eq.  D.l) 

yields  for  p  ~  exp  (-iuvt) 

V2P  +  k(z)2p  =  0  (Eq.  D.  2) 

where  k(z)  is  the  wavenumber. 

The  "eikonal-solution" 

p  =  A(r,z)  exp[ i  S(r,z)]  (Eq.  D.3) 


permits  for  slowly  varying  A(r,z)  the  "eikonal  equation" 


(VS)2  =  k(z)2 


(Eq.  D. 4) 


For  the  derived  solution  for  p  the  phase  function  S(r,z)  permits  the 
surfaces  of  constant  phase,  and 


S(r,z)  =  const.  (Eq.  D.5) 


are  the  wavefronts  of  the  sound  field.  The  vector  VS  are  normal  to  the 
wavefront  and  chosen  to  form  tangents  to  space  curves,  the  so-called 
"rays". 
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The  path  element  along  a  ray  in  range  and  depth  coordinates  (r,z)  is 
given  by: 

d £  =  J  1  +  (dr/dz)2'  dz  (Eq.  D.6) 

and  the  angle  0  which  the  ray  makes  with  the  horizontal  can  be  found  as 

^|  =  sin  6  ;  =  cos  0  ;  ^  =  cot  0  (Eq.  D.7) 


The  unit  vector  pointing  normal  to  the  surface  S(r,z)  =  const,  and 
being  parallel  to  the  vector  VS  is  defined 


n  =  (dr/d£,  dz/d£) 

Together  with  Eq.  D.4  follows 

dr  _  8S  ,  dz  _  §S 
K  d£  9r  *  K  d£  3z 


with  the  operator  (d/d £) 

d_  /,  dr  \  _  a  dS 
d£  y  d£  J  3r  d £ 

and  with 

dS  as  dr  t  3S  dz 
d£  3r  d£  3z  d£ 

one  has  by  using  Eq.  D.9 


with  Eq.  D.10  again 

3l(k5l)=s(k  cos  e)  =  0 

Going  along  d SL  ,  i.e.  a  ray,  one  has 

k(z)  cos  0  =  const,  or  the  Snell1 

cos  0/c(z)  =  cos  0  /c  =  const. 

oo 


(Eq.  D. 8) 

(Eq.  D.9) 

(Eq.  D. 10) 

(Eq.  D. 11) 

(Eq.  D. 12) 

law  (Eq.  D. 13) 

(Eq.  D. 14) 


This  law  determines  the  variation  of  grazing  angle  0  (the  angle  relative 
to  the  horizontal)  along  the  ray  corresponding  to  the  variation  of  sound 
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speed  with  depth  z  ,  and  relating  it  to  the  corresponding  quantities  0Q, 

C  at  the  source.  From  this  and  Eqs.  D.6,  D.7  the  ray  equation  is  derived 
o 

r(z)  =  f  dZ . -  —  7  (Eq.  D.  15) 

zq  Vco/c(z)  •  cos  0)2-l 

for  a  ray  coming  from  the  source  at  a  grazing  angle  0Q  .  Solving  this 
equation  for  a  constant-gradient  profile,  yields: 


r(z)  =  f  -~~=  -  Ja2  -  z2'  (Eq.  D.16) 

z  V«2-z2 
with  ° 


c(z)  =  bz 


and 

cos  0 

o 


This  permits 

a2  =  r2  +  z2  (Eq.  D.17) 


with  proper  choice  of  integration  constants 

Hence  for  a  small  gradient,  the  "rays"  are  circles  with  large  radius,  for  a 
large  gradient  they  are  circles  with  small  radius  and  for  isovelocity 
profile  the  rays  become  straight  lines.  This  concept  can  explain  the  sound 
channeling  around  the  SOFAR  channel  axis. 

Effects  of  boundary  occur,  when  upward  curving  rays  are  reflected  by  the 
sea  surface  and  downward  rays  are  reflected  by  the  sea  bottom.  Bottom 
reflections  are  treated  by  this  concept  assuming  equal  grazing  angles  of 
incidence  and  reflections,  but  with  reduced  amplitudes  after  reflection. 
The  sea  surface  is  almost  a  perfect  reflector,  for  a  low  frequency 
(<  100  Hz),  but  the  phase  is  changed  by  n  after  the  reflection.  If  a  ray 
has  touched  a  caustic,  which  is  a  curve  forming  a  geometrical  envelope  of 
rays,  the  phase  is  changed  by  n/2.  In  fact  modified  ray  theory  does 
consider  the  caustic  correction. 


In  calculating  the  phase  explicitly  between  source  and  receiver  along  the 
rays,  Eq.  D.4  yields  by  integrating  over  a  ray-path. 


S  =  /  kd  £=  u)J^  =  uj/dt  =  uoT 


(Eq.  D. 18) 


where  T  is  the  travel  time  along  the  ray. 


Each  ray  amplitude  has  to  have 
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a  phase  factor  exp(iuiT)  for  coherent  summation.  In  writing  Eq.  D.18 
expl icitly 


u»T  =  u>  J 


dz 


c  sin  6 


=  «*>  S 


dz 


c  sin  0 


(sin  02  +  cos2  0) 


f  sin  0  ,  .  r  cos2  0  , 

=  u)  f  -  dz  +  lu  J  - 7 — „dz 

J  c  J  c  sin  0 


(Eq.  D. 19) 


Ix  =  ui  f  ~~  /  1  ■  cos2  0 


(Eq.  D. 19a) 


and  with 


cos  0  =  n  cos  0 
o 

n  =  co/c(z) 


k  =  u)/c 
o  o 


we  have 


I,  =  k  J  7  n2  -  cos2  0  dz 

a  n  J  *  n 


(Eq.  D. 19b) 


For  the  second  integral 

dz 


I2  =  «*»  J 


c  si n  0 


cosz0 


2fi  l 


(Eq.  D. 20a) 


follows  with,  Eq.  D.7,  dz/d£=  sin  0 
T  _  c  cos2  0 

I2  -  “J - — 


(Eq.  D. 20b) 


With  the  horizontal  component  k^  =  k  cos  0  of  the  wavenumber  k  =  w/c  one 
finds  from  Snell's  law  the  horizontal  component 

k  cos0  =  k  cos  0  =  const, 

o  o 

being  range- independent,  i.e.  a  constant  along  the  ray.  This  yields  for 
Eq.  D. 20b) 

I3  =  J  k  dr  =  (k  cos  0  )r  (Eq.  D.20c) 

J  r  0  0 

The  total  phase  along  the  ray  can  be  expressed 

ujT  =  (k  cos  0  )  r  +  k  f  J  n2-cos2  0  dz 
0  0  0  J  y 

^0 

=  J  krdr  +  /  kzdz  (Eq.  D.21) 
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where  the  vertical  wavenumber 

k  =  k  sine  =  —  /n2-cos2  6 
z  c  ^ 

o 

is  depth  dependent  through  n  =  cq/c(z) 

From  this  the  well  known  relation 

k2  =  k2+  kz  (Eq.  D. 22) 

is  found. 


D. 2  NORMAL  MODE 

With  a  monopole  source  at  r  =  0  and  z  =  z  the  normal  mode  soluti 
for  the  wave  equation  is  0 

6(r)  6(z-z  ) 


on 


V2p  +  kp  =  - 


2  nr 


(Eq.  D. 23) 


with  separation  of  variables,  the  solution  is  of  the  form 


p(r,z)  =  I  Rn(r)  Dn(z) 
n 


(Eq.  D. 24) 


with  some  constants  a^  and  k^  we  can  take  the  range  function 
Rn<p)  =  an14Ho1)  (knr> 


(Eq.  D. 25) 


which  is  the  Green's  function  of  the  wave  equation,  satisfyi 


ng 


11-  (r-^V  k'2  R  -  -  «(r) 
r  dr  y  dr  J  K  n  Kn  2ztr  n 


and  for  cylindrical  symmetry  one  has 

V2p  =  I  <L  L  §|\+  sie 

r9rl  r9/  Oz2 


(Eq.  D. 26) 


(Eq.  D. 27) 
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Inserting  Eqs.  D.25  and  D.26  in  Eqs.  D.23  and  D.24  yields 


r  d2D  2  1 

I  -k2  R  D  +  R  -r£  +  ^  R  D 
I  n  n  n  n  dz2  c2  n  nj 


<5(r)  6(z'zn) 

=  Ml)  z  a  D  (Z) - ^ - 2 

„  n  n  2nr 

2nr  n 

The  right  hand  side  vanishes  for 

a  =  D  (z  ) 
n  n  o 

and  using  for  depth-function  Dn(z) 

I  D  (z)  D  (z  )  =  6(z-z  ) 
n  n  o  o 

n 


This  function  satisfies  the  depth  equation. 


(Eq.  D. 28) 

(Eq.  D. 29) 

(Eq.  D. 30) 


d2  °n 

dz2 


-  k2  )  D 
n  I  n 


=  0 


(Eq.  D. 31) 


This  "eigenvalue  equation"  permits  solutions  for  its  characteristic 
"eigenvalues"  kn  under  the  constraint  that 


D  (z  )  =  D  (z  )  =  0 
n  s  n  oo 


(Eq.  D.  32) 


The  eigenvalue  kn  forms  a  "discrete  spectrum"  or  a  continuum. 
The  solution  for  an  outgoing  wave  is  now  obtained  from  Eq.D.24. 


.  .  oo  (i) 

P(r.z)  =  e‘lujt  l  I  H  (k  r)  Dn(z  )  0  (z)  , 

n=l 


(Eq.  D. 33) 


known  as  normal  mode  series.  At  large  ranges  r  » 


H(1)  (k  r) 
o  v  n 


/ 


"  Kk  r  -  J) 


7t  knr 


(Eq.  D.  34) 


corresponds  to  cylindrical  waves, 

showing  that  the  horizontal  wave-numbers  are  constant,  as  they  are  in  ray 
theory,  too.  In  other  words,  they  are  range-independent.  The  ray-mode 
equivalence  in  terms  of  an  "equivalent  ray  angle"  may  help  interpretations . 
The  horizontal  wave  number  kp  for  the  mode  solution  is  assumed  to  be 

equal  to  the  horizontal  wave  number  kr  =  w  cos  6/c  in  ray  theory. 
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t  h 

Hence  the  n  normal  mode  might  be  associated  with  an  "equivalent  ray  angle 
6"  of  that  mode  by 


.1  k  c(z) 

en(z)  =  cos  — -  ,  (Eq.  D.35) 

showing  the  difference  that  at  given  depth  z  ,  arrivals  in  terms  of  rays 
are  continuous  and  in  terms  of  modes  are  a  discrete  set  of  arrivals. 
Introducing  a  vertical  wavenumber  one  obtains  a  depth  equation, 

Dn(z)  +  kj  Dn(z)  =  0  (Eq.  D.36) 

which  can  be  solved  analytically  for  certain  profiles. 

The  Eq.  D.31  might  be  rewritten  by  introducing  the  "potential  function" 


V(z)  =  k2  -  k2(z) 


(Eq.  D. 37) 


with  kQ  =  w/co  at  channel  axis  and  the  new  eigenvalue 
E  =  k2  -  k2 

non  (Eq.  D. 38) 


the  form  of  the  depth  equation  is  found 
d2  Dn(z)  r  1 

~dP -  +  LEn  "  V(Z)J  Dn(z)  =  0  (Eq.  D-39) 

Hence  for  a  given  sound  speed  profile  a  related  potential  function  exists. 
The  infinity  value  of  this  potential  determines  the  boundary  between 
continuous  and  discrete  mode  propagation  and  is  given  by 


(Eq.  D. 40) 


The  turning  point  for  mode  functions  occurs  if  —  V(z)  ,  whereas  beyond 

this  depth  z  the  vertical  wavenumber  is  negative  and  the  depth  function 
becomes  exponentially  decaying.  For  E  >  V  the  waves  are  always 

travelling  towards  z  ->  »  ,  hence  leaking  into  the  bottom.  This  is  known 
as  "leaky  modes". 

The  mode  cut-off  occurs  when  Ep  >  iu2/c2  ,  which  is  higher  than  the  margin 

of  the  continuum  since  the  latter  is  qiven  bv  E  >  V 

M  J  n  a>  • 

From  this  short  description  it  is  generally  concluded  that  in  shallow  water 
the  boundary  conditions  at  the  bottom  might  cause  an  early  cutoff,  and  only 
few  modes  propagate.  In  deep  water  at  very  low  frequencies  only  few  modes 
are  trapped  in  the  "potential  pot".  In  both  cases  ray  theory  cannot 
describe  accurately  the  acoustic  propagation. 
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APPENDIX  E 

LISTINGS  OF  THE  COMPUTER  PROGRAMS  USED  FOR  THE  SIGNAL  PROCESSING 


E. 1  MAIN  PROGRAM  FOR  OFF-LINE  PROCESSING 

The  purpose  of  the  off-line  processing  is  described  in  Chapter  3.  The 
following  is  a  listing  of  its  main  program. 


0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 

0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 


FTN4, L 
*EMA(RBE,  5) 

PROGRAM  LUNDE (3, 99) ,  ****  SPG  BEAMFORMING  3/3/81  ****** 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  PROGRAM  FOR  CONVENTIONAL  AND  MAXIMUM  LIKELIHOOD  BEAMFQRMING  C 
C  COMBINED  WITH  DIFFERENT  METHODS  FOR  OWN  VESSEL  NOISE  REDUCTION  C 
C  C 
C  DVELOPED  BY  E,  B.  LUNDE  (SPG)  C 
C  MARCH  1981  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  ARRAYS  ARE  DEFINED  SO  THAT  THERE  CAN  BE  MAXIMUM 

C  19  ARRAY  HYDROPHONES  +  1  REFERENCE  HYDROPHONE, 

C  181  BEAMS,  10  FREQUENCIES,  AND  10  PROCESSING  MODES 

C 

C  PARAMETER  VALUES  : 

C 

C  NH  -  NUMBER  OF  HYDROPHONES 

C  NB  -  NUMBER  OF  BEAMS 

C  NF  -  NUMBER  OF  FREQUENCY  CELLS 

C  F  -  ARRAY  WITH  FREQUENCY  VALUES 

C  C  -  SOUND  SPEED  (m/sec) 

C  A  -  HYDROPHONE  SEPARATION  (m) 

C  NA  -  NUMBER  OF  AVERAGE  IN  THE  RAW  MATRIX 

C  NT  -  NUMBER  OF  TOTAL  AVERAGES  PER  TAPE 

C 

C  NH, NB, NF, F, A,  C  PARAMETERS  WILL  BE  SET  BY  THE 
C  ACQUSITION  PROGRAM  AND  WILL  BE  IN  THE  MAGNETIC 

C  TAPE  HEADER  BLOCK. 

C 

C  PARAMETERS  NA  ,  NT  WILL  BE  TYPED  BY  THE  OPERATOR  . 

C 

C  V.  DUARTE  (COM)  05/03/81/ 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

DIMENSION  R  <  400, 2 ) , E ( 400, 2), EV(20,  4) 

DIMENSION  BGR ( 182 ) , PR ( 20,  6 ) 

DIMENSION  B ( 20, 2), BE(20, 2), BSQC20) 

DIMENSION  NN<5>, IDCBC144) 

DIMENSION  TAB (181,2) 

DIMENSION  P(182, 10), PSP(6) 

DIMENSION  F< 10) , KM AX ( 20 ) ,  H ( 20,  2,  2) 

DIMENSION  IB(815), FIB<12),  RB(1 83 ) 

DIMENSION  RSP (36, 2), BSP (11,  11) 

EQUIVALENCE  (P  ( 1,  1 ) ,  E<  1,  1 ) ,  IOCBCD).  (P  ( 1,  6) ,  R  ( 1,  1  )  ) 

EQUIVALENCE  ( IB ( 1 ) , RB < 1 ) , H( 1 , 1 ,  1 ) ) ,  ( IB ( 19 ) , FIB ( 1 ) ) 

EQUIVALENCE  ( IB ( 161 ) , B < 1 ,  1 ) ) ,  ( IB ( 241 ) , BE ( 1 ,  1)) 

EQUIVALENCE  ( IB (321 ), BSQ(1 ) >, ( IB < 361 ) , EV ( 1 , 1 ) ) 

EQUIVALENCE  ( IB ( 521 > , BGR ( 1 ) ) , (BSP<1, 1),PR<1,1>) 

C 

COMMON/RBE/LREC (15),  TREC (400) , RF<2, 4000), 

*BEW ( 4000, 2 ) ,  BSQW ( 4000 ) 
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0056 

0057 

0058 

0059 

0060 

0061 

0062 

0063 

0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 

0081 

0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 

0100 

0101 

0102 

0103 

0104 

0105 

0106 

0107 

0108 

0109 

0110 

01 1 1 

0112 

01 13 

0114 

0115 

0116 

0117 

0118 

0119 

0120 

0121 

0122 

0123 

0124 
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c 

DATA  NP/10/ 

C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c 

c 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


10  CONTINUE 

CALL  EXEC (3, 407B > 

CALL  EXEC  < 1 1 7,  IB, 527) 

WR ITE< 1 1  1000)  <  IB ( I  ) »  1  =  10,  13) 

1000  FORMAT ( “  TAPE  HEADER  :  ", 4A2, /, 

READ ( 1 , 1001 )  IP 

1001  FORMAT (2A2) 

IF ( IP. EG. 2HYE )  GO  TO  30 
WRITEd,  1002) 

1002  FORMAT ( /, "  INSERT  NEW  TAPE  AND  TYPE 
20  CONTINUE 

READ< 1,1001)  IP 
IFdP.EQ  2HG0)  GO  TO  10 
GO  TO  20 
30  CONTINUE 

WRITEd,  1035) 

1035  FORMAT ( /,  "  MAG  TAPE  NUMBER  ?  _"  ) 
READ ( 1 , * )  MGTP 
NH=IB ( 16 ) -1 
NB=IB ( 17 ) 

NF-IB(18) 

A=FIB(NF+1 ) 

C=F IB ( NF+2 ) 

WRITEd,  1003)  NH,  NB,  NF,  A,  C 

1003  FORMAT </,"  NUMBER  OF  HYDROPHONES  = 


OK  ?  < YE/NO) 


"GO"  TO  PROCEED 


12,  /, 


*"  NUMBER  OF  BEAMS  =  ",  13, /, 

*"  NUMBER  OF  FREQUENCIES  »  ",  12,  /, 

$"  HYDROPHONE  SEPARATION  =  ",F6.  4,  /, 

*"  SPEED  OF  SOUND  (m/s>  -  ", F12. 4, /> 

DO  40  1=1, NF 
F ( I ) =FIB  < I ) 

WRITEd,  1004)  I,  F<  I  > 

1004  FORMAT (2X,  12,  4X,  F8.  4,  "  Hz  ") 

40  CONTINUE 

WRITEd,  1005) 

1005  FORMAT ( /, "  OK  ?  (YE/NO)  _"  > 

READ<1, 1001)  IP 

IFUP.NE  2HYE)  STOP  001 
4141  CONTINUE 

WRITE< 1, 1006) 

1006  FORMAT ( /  ,  11  ENTER  STARTING  BLOCK  NUMBER  ON  M.  T. 
READ ( 1 , *  >  IC 

IF( IC  GT. 11400)  GO  TO  4141 
1106  CONTINUE 

WRITEd,  1007) 

1007  FORMAT ( /,  "  NUMBER  OF  AVERAGES  ?  _•' > 


READ  d , * )  NA 

IF ( NA.  GE.  20)  GO  TO  50 

WRITEd,  1107) 

1107  FORMAT ( /,  "  INPUT  ERROR  /,  “  20<=  NA  <=600  "> 

GO  TO  1106 

50  CONTINUE 

WRITEd,  1008) 

1008  FORMAT (/,"  NUMBER  OF  AVERAGES  PER  TAPE  ?  _" > 

READ  d ,  * )  NT 

IF(NT,  LE.30.  AND  NT  GT.O)  GO  TO  1108 
WRITEd,  1118) 

1118  FORMAT < /,  "  INPUT  ERROR  0<  NA*NT  <=  600  " ) 

GO  TO  1106 

1108  CONTINUE 


) 
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0125 

0126 

0127 

0128 

0129 

0130 

013.1 

0132 

0133 

0134 

0135 

0136 

0137 

0138 

0139 

0140 

0141 

0142 

0143 

0144 

0145 

0146 

0147 

0148 

0149 

0150 

0151 

0152 

0153 

0154 

0155 

0156 

0157 

0158 

0159 

0160 

0161 

0162 

0163 

0164 

0165 

0166 

0167 

0168 

0169 

0170 

0171 

0172 

0173 

0174 

0175 

0176 

0177 

0178 

0179 

0180 

0181 

0182 

0183 

0184 

0185 

0186 

0187 

0188 

0189 

0190 

0191 

0192 

0193 

0194 

0195 

0196 

0197 
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c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C  WRITE  ON  FILE  "  >LUNDE  11  THE  FOLLOWING  INFO.  : 

C 

C  NT  -  NUMBER  OF  AVARAGES  PER  TAPE 

C  NF  -  OF  FREQUENCIES 

C  NP  —  11  OF  PROCESSINGS 

C  FI  -  FREQUENCY  #  1 

C  :  :  :  : 

C  Fn  -  H  n  FOR  n=NF 

C 

C  THIS  FILE  WILL  BE  READ  BY  THE  PLOTTING  PROGRAM 

C 

C 

CALL  QPEN< IDCB,  IERR,  6H>LUNDE) 

IF ( IERR.  LT.  0 )  STOP  100 
CALL  CODE 

WRITEdB.  1109)  NT i  NF/  NP 
1 109  FORMAT ( 2Xi  314) 

CALL  WRITFdDCB,  IERR,  IB,  8) 

IF < IERR.  LT. 0)  STOP  101 
CALL  CODE 

WRITEdB,  1119)  (F<I),  1  =  1,  10) 

1119  FORMAT ( 10(2X,  F6.  2) ) 

CALL  WRITFdDCB,  IERR,  IB,  NF*4) 

IF ( IERR.  LT. O )  STOP  102 
CALL  CLOSEdDCB) 

C 

C 

C  END  OF  INPUT  PARAMETER  PHASE  CVD1. 

C 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  SOME  USEFUL  CONSTANTS 

C 

NFE=NF+1 

NBE*NB+1 

NBR=NB~1 

NBH=NB/2+l 

NHE=NH+1 

NHD=NH*2 

NHED=NHE*2 

NHSQ=NH*NH 

NHESG=NHE*NHE 

NSP=NH/2 

ALN1G=ALGG( 10. ) 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C  CALCULATE  TABLE  USED  FOR  STEERING  VECTOR  GENERATION 

C 

P I  =4.  *ATAN<  1.  ) 

ARG=2.  *PI/U.  *NBR> 

SNH=1 .  /SORT ( 1 .  *NH) 

DO  81  JB  =  1  #  NB 
JBM=UB-1 

TABCUB, 1 >=COS(ARG*JBM)*SNH 
TAB ( JB i 2 ) =S I N ( ARG*UB  M  > *SNH 

81  CONTINUE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  FIND  THE  MAG  TAPE  START  POINT 

C 
C 

IF(IC.  LE.  1 )  GO  TO  83 

82  CONTINUE 

CALL  EXEC  C 1 ,  7,  IB, 815) 

IF ( IB ( 1 ) .  EQ.  IC)  GO  TO  83 
GO  TO  82 

83  CONTINUE 
C 
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0198 

0199 

0200 

0201 

0202 

0203 

0204 

0205 

0206 

0207 

0208 

0209 

0210 

0211 

0212 

0213 

0214 

0215 

0216 

0217 

0218 

0219 

0220 

0221 

0222 

0223 

0224 

0225 

0226 

0227 

0228 

0229 

0230 

0231 

0232 

0233 

0234 

0235 

0236 

0237 

0238 

0239 

0240 

0241 

0242 

0243 

0244 

0245 

0246 

0247 

0248 

0249 

0250 

0251 

0252 

0253 

0254 

0255 

0256 

0257 

0258 

0259 

0260 

0261 

0262 

0263 

0264 

0265 

0266 

0267 
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cccccccccccccccccccccccccccccccc 

c 

C  START  OF  MAIN  LOOP 

C 

DO  191  JT»=1  /  NT 
C 

cccccccccccccccccccccccccccccccccccccccccccccc 

c 

C  CLEAR  THE  CORRELATION  MATRICES  IN  EMA 

C 

CALL  WSMY  <  0.  ,  RF ( 1 ,  1 ) ,  1 , RF  < 1  *  1 > ,  1 , 8000  > 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CCC 

CCC  READ  RECORDS  FROM  TAPE 
CCC 

DO  101  JA=1 ,  NA 
CALL  LEMA1 (7, I ERR ) 

IF ( I ERR  LT.  0)  STOP  555 

CCC 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

c  SUBTRACT  1/48  OF  X(NH)  FROM  X(NHE),  THIS  BECAUSE  OF 

C  THE  METHOD  USED  IN  THE  BEAMFORMER  TO  DELAY  X  <  NHE ) 

C 

G— 1.  /48. 

DO  22  KA=1,2 

LM=NHD+KA 

LH=LM-2 

CALL  WP I V ( G,  TREC ( LH ) ,  NHEDi  TREC (LM), NHED,  TREC  <  LM ) ,  NHEDi  NF ) 

22  CONTINUE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  CALCULATE  THE  OUTER  PRODUCT  OF  THE  VECTOR  X 

C  AND  ADD  THE  RESULTING  MATRIX  TO  MATRIX  RF 

C 

DO  23  JF=1, NF 
LF= ( JF- 1 >*NHED+1 
MF= ( JF- 1 ) *NHESQ+ 1 

CALL  WVMOV <  TREC ( LF )  ,  1,R<1,  1),  1 , NHED ) 

DO  24  JH=1,NHE 
LH= ( JH-1 ) *2 
LFH=LF+LH 
MFH=MF+<UH-1 )*NHE 
YR=R ( LH+ 1 , 1 ) 

YI=R<LH+2, 1) 

CALL  WPIV<YR,  TREC<LF>,  1,RF<1,MFH),  1,RF(1,MFH),  l.NHED) 

CALL  WP I V<  YI i  TREC  <  LF+ 1 ) ,  2,RF(1, MFH ) #  2/ RF<1,  MFH ) *  2» NHE) 

CALL  WPIVC-YI, TREC(LF),  2,  RF(2,  MFH),  2,  RF(2,  MFH),  2,  NHE) 

24  CONTINUE 

23  CONTINUE 
101  CONTINUE 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  START  OF  FREQUENCY  LOOP 

c 

DO  21  JF=1 ,  NF 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  READ  THE  CORR  MATRIX  AT  THE  ACTUAL  FREQUENCY  FROM  EMA 

C  TO  THE  MATRIX  R,  NORMALIZE  AND  MOVE  IT  TO  MATRIX  E 

C 

LF=(UF-1 )*NHESQ+1 
DO  115  KA«1,  2 

CALL  WVMOV ( RF < KA,  LF),  2,  R(l,  KA>,  1, NHESQ) 

115  CONTINUE 

CALL  VSUM ( TR , R  < 1 ,  1 ) , NHE+ 1 , NH ) 
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0268 

0269 

0270 

0271 

0272 

0273 

0274 

0275 

0276 

0277 

0278 

0279 

0280 

0281 

0282 

0283 

0284 

0285 

0286 

0287 

0288 

0289 

0290 

0291 

0292 

0293 

0294 

0295 

0296 

0297 

0298 

0299 

0300 

0301 

0302 

0303 

0304 

0305 

0306 

0307 

0308 

0309 

0310 

0311 

0312 

0313 

0314 

0315 

0316 

0317 

0318 

0319 

0320 

0321 

0322 

0323 

0324 

0325 

0326 

0327 

0328 

0329 

0330 

0331 

0332 

0333 

0334 

0335 

0336 

0337 

0338 

0339 
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P0W=TR/(i. *NH) 

P0WI-1. /pow 
DO  111  JH*1#NHE 
LH=( JH-1 >*NHE+1 
MH=( JH-1 )*NH+1 

CALL  VSMY(POWI, R<LH*  1 >.  1j  E(MH»  1 >.  1»  NH> 

CALL  VSMY (-POWI , R(LH»  2) »  1.  E<MH»  2),  i,  NH) 

113  CONTINUE 
111  CONTINUE 

POW=PGW/< 1.  *NA) 

QS=R(NHESQ,  1>*P0WI 

CALL  VSMY  ( 100.  ,  E<1,  1  ),  NHE,  PR  ( 1 , 2 ) ,  1 ,  NH> 

PR ( NHE* 2)=QS*100. 

PR ( NHE, 1 )=POW 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C  PRELIMINARY  SPLIT  BEAM  PROCESSING 

C 

NN( 1 )=-NH+l 

NN ( 2 ) =“NH+NSP*NH+NSP+1 
NN  ( 3  >  *=-NH+NSP  #NH“NSP-1 
NN  ( 4 ) =-NH+NSP *NH+ 1 
DO  701  JSP35!  i  4 
LINC=NHE 

IF  (JSP.EQ.3)  LINC^-NHE 
LSP*  <  JSP-1 >*NSP 
DO  702  JH=1 / NSP 
NN ( JSP  >  =NN ( JSP ) +NH 
NEL*»NSP  +  1-JH 
IF  (JSP.EQ.3)  NEL-JH— 1 
DO  703  KA=1,2 

CALL  VSUM ( RSP ( JH+LSP/  KA ) ,  E ( NN ( JSP ) ,  KA ) * LINC, NEL) 

703  CONTINUE 
702  CONTINUE 
701  CONTINUE 

C  WR ITE( 6/ 822)  ( (RSP (JSP/ KA), KA«1, 2), JSP=1, 36) 

822  FORMAT (2F10. 3) 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  USE  A  SUBROUTINE  TO  FIND  EIGENVALUES  AND  EIGENVECTORS  OF  R 

C  THE  EIGENVALUES  ARE  STORED  IN  DECREASING  ORDER  IN  VECTOR  EV 

C  AND  THE  RESPECTIVEIGENVECTQRS  IN  MATRIX  E 

C 

NH1=NH 

CALL  EIGEN  (E<  1 ,  1),E(1,2),R(1,  1  ) ,  R  ( 1 ,  2 ) ,  NHSQ,  EV  ( 1 ,  1),NH1) 

CALL  VSMY (10.  , EV< 1 ,  1 ) ,  1, PR ( 1.  1 ) ,  1#  NH> 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C  INVERT  THE  EIGENVALUES 

C 

CALL  VSDV ( 1 .  ,  EV(1#  1>«  1iEV(1j2>»  ItNH) 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C  START  OF  PRELIMINARY  BEAM  LOOP 

C 

DO  41  JB=1 #  NBE 
IF  (JB.  EQ.  NBE)  GO  TO  223 
C 
C 

C  CREATE  STEERING  VECTOR  FROM  TABLE 

C 

NBH=NB/2+l 

K=NBH-JB 

IF  (K.  LT.  0)  K-K+NBR 
L=0 

DO  33  JH=1#  NH 

DO  32  KA=1,2 

B( JH,  KA )=TAB (L+l/  KA) 


65 


0340 

0341 

0342 

0343 

0344 

0345 

0346 

0347 

0348 

0349 

0350 

0351 

0352 

0353 

0354 

0355 

0356 

0357 

0358 

0359 

0360 

0361 

0362 

0363 

0364 

0365 

0366 

0367 

0368 

0369 

0370 

0371 

0372 

0373 

0374 

0375 

0376 

0377 

0378 

0379 

0380 

0381 

0382 

0383 

0384 

0385 

0386 

0387 

0388 

0389 

0390 

0391 

0392 

0393 

0394 

0395 

0396 

0397 

0398 

0399 

0400 

0401 

0402 

0403 

0404 

0405 

0406 

0407 

0408 

0409 

0410 

0411 

0412 


SACLANTCEN  SM-158 


32  CONTINUE 
L=L+K 

IF  <L. GE. NBR)  L=L-NBR 

33  CONTINUE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  CALCULATE  BEAM  DIRECTION  OR  RELATIVE  WAVENUMBER 

C 

AL=A*F(JF)/C 

JJ=1 

ARG= ( NBH-JB ) / ( AL*NBR ) 

225  CONTINUE 

AARO=ABS ( ARO ) 

IF  (AARG.  LE.  1.  )  GO  TO  222 
BGR ( JB ) =-10. *AARG 
GO  TO  223 

222  CONTINUE 
ARH-SGRT ( 1 .  -ARG*ARG) 

BGR  ( JB  )  =:ATAN2<  ARH,  ARG)*180.  /PI 

223  CONTINUE 

C  IF  <JJ.  EG.  2)  GO  TO  224 

C  JJ=2 

C  IF  (JB.LE.  NBH)  ARG=ARG-i.  /AL 

C  IF  (JB.GT.  NBH)  ARG=ARG+1.  /AL 

C  GO  TO  225 

C  224  CONTINUE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C  CALCULATE  THE  VECTORS  B'E,  B 'G,  G'E,  AND  G'G,  AND  STORE  IN  EMA 

C 

IF  (JB.NE.  NBE)  GO  TO  47 
DO  48  KA«1,2 

CALL  VMOV ( E < NHSG+ 1 ,  KA),  1, B< 1, KA) ,  1,  NH) 

C  NSB  (NBE,  KA)*0 

48  CONTINUE 
47  CONTINUE 

LB=(JB-1)*NHE+1 
DO  42  JH-1,  NHE 
LH= ( JH— 1 ) *NH+ 1 
DO  43  KA=1#  2 
DO  44  KB-l ,  2 

CALL  VDOT ( H ( JH/  KA#  KB ) #  B  < 1 #  KA ) ,  1,  E<LH>  KB ) ,  1 »  NH) 

44  CONTINUE 
43  CONTINUE 
42  CONTINUE 

CALL  VADD ( H  ( 1 ,  1 ,  1 ) ,  1 , H( 1 ,  2,  2 ) z  1 , BE( 1/  1 )  i  1, NHE) 

CALL  VSUB  (H(  1  j  2.,  1 )  #  1*  H(  1«  1#  2)#  If  BE(  1#  2)  i  1,  NHE) 

DO  45  KA=1  /  2 

CALL  VWMOV ( BE ( 1 # KA),  1, BEW(LB,  KA),  1, NHE) 

45  CONTINUE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  CALCULATE  THE  SQUARE  AMPLITUDE  OF  EACH  ELEMENT  OF  MATRIX  BE 

C  AND  STORE  IN  EMA 

C 

DO  46  KA=1 , 2 

CALL  VMPY  ( BE  <  1 ,  KA),  1,  BE<  1,  KA),  1,  H(l,  KA,  1  ),  1 ,  NHE ) 

46  CONTINUE 

CALL  VADD(H( 1,  1 ,  1 ) ,  1 , H ( 1 , 2,  1 ) ,  1 , BSG( 1 ),  1, NHE) 

CALL  VWMOV ( BSG ( 1 ) , 1, BSGW(LB), 1, NHE) 

P ( JB, 7 ) =100.  *BSQ(i) 

P( JB, 8 )  =  1 00.  *BSQ(2) 

P< JB, 9 )  =  100.  *BSGCNHE> 

IF  (JB.LT.  NBE)  GO  TO  715 
GNRM*SQRT  <  BSG  C  NHE ) ) 

QNRMI  =  1.  /GNRM 

CALL  VSMY ( GNRM I , P ( 1 , 9 ) ,  1 , P ( 1 , 9 ) ,  1 , NB ) 

CALL  VSMY ( GNRM  1*100.  ,BSG(1),  1,PR(1,  5),  1,  NHE) 

715  CONTINUE 

C  WRITE<6, 819)  JB, WA, WB,  <BSQ( JH) , JH=1, NHE) , JB 
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C  819  FORMATCX,  13,  2F8.  1,  2X,  20F4.  3,  2X#  16) 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C  CONTINUE  SPLIT  BEAM  PROCESSING 

C 

DO  711  JSP=1,  3 
NEL=NSP-1 

IF  CJSP.EQ.  3)  NEL=NSP*2-1 
LSP= ( JSP-1 )*NSP+1 
DO  712  KA=1, 2 
DO  713  KB=1 , 2 

CALL  VDOTCHC JSP, KA,  KB),  BC2,  KA),  1, RSP ( LSP+1 , KB ) ,  1, NEL ) 

713  CONTINUE 
712  CONTINUE 

PSPC JSP)=RSPCLSP,  1  )/( 1.  *NH)+2.  *CHC JSP,  1,  1 )-HC JSP,  2,  2) ) *SNH 
711  CONTINUE 

PSP  C3)=PSP  C3) /2. 

PSP  C  4 )  =  CH  C3, 2i  1 )+H(3i  1«  2) >*SNH 
PSP ( 5) =PSP  <3) /SORT  (PSP  C 1 )*PSP (2) ) 

PSP ( 6 ) »PSP ( 1 ) +PSP ( 2 >  +2.  *PSP(3) 

P< JB, 6)=PSP(3)*100. 

PC  JB,  10)=1. 

IF  (BGRC  JB).  LT.  0.  )  GO  TO  718 
DEL=ATAN2 ( PSP  C  4  > ,  PSP (3) >/<2.  *PI*AL*NSP> 

DCOS-COS ( BGR 80.  ) 

RCOS=DCOS-DEL 

IF  ( ABS(RCOS) .  GT.  1 .  )  GO  TO  718 
RSIN=SGRT ( 1 .  -RCOS*RCOS> 

P  <  JB,  10)=ATAN2<RSIN, RC0S)*180.  /PI 

718  CONTINUE 

IF  ( ISSW<3).  GE.  0)  GO  TO  719 

WRITE <6, 817)  JB# BGR ( JB > *  <PSP< JM).  JM=1# 6 ) , JB 
817  FORMAT <  X,  13,  4X,  F10.  2,  4X,  6F10  4,  4X#  16) 

719  CONTINUE 
C 

41  CONTINUE 
C 

C  END  OF  PRELIMINARY  BEAM  LOOP 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  WRITEC6, 818)  JA# JF, POW, QNRM, QS 

C  818  FORMAT <2110, 3F10. 3 ) 

IF  ( ISSW<3).  LT.  0)  WRITEC6, 535) 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  FIND  THE  BEAM-INDEX  ASSOCIATED  WITH  EACH  EIGENVALUE 

C 

DO  51  JH=1,  NHE 

CALL  WMAX ( KMAX  <  JH ) , BSQW < JH ) , NHE, NB ) 

PRC JH,  3 ) =KMAX ( JH ) 

PR ( JH,  4 ) =BGR ( KMAX ( JH ) ) 

51  CONTINUE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  PRINT  THE  EIGENVALUES  AND  THE  ASSOCIATED  BEAM-INDECES 

C 

C  WRITEC6,  810)  <JH,JH=1,NH) 

WRITEC 1,810)  ( JH, JH=1 , NH) 

C  WRITEC6, 810)  ( KMAX ( JH ) , JH=i , NH ) 

WRITEC 1, 810)  CKMAXC JH), JH«1# NH) 

C  WR ITEC6, 810 )  C PR C JH, 1 ), JH=1 # NH ) 

WRITEC 1, 810)  CPRC JH, 1 ) , JH=1 , NH ) 

810  FORMATCX, 1914) 

811  FORMATCX,  F5.  3,  18F4.  3) 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C  FIND  REF  BEAM  DIRECTION 
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BDEL=-1 . 

DO  56  JB=1,NB 

IF  (BGR< JB).  LE.  0.  )  GO  TO  56 
ADEL=P ( JBj  lO)-BGR(JB) 

IF  (ADEL.  LE.  0.  .  AND.  BDEL.  GT.  0.  )  GO  TO  57 
BDEL=ADEL 

56  CONTINUE 

57  CONTINUE 
KREF=JB 

IF  (BDEL. LT.  -ADEL)  KREF=JB-1 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C  FIND  THE  EIGENVALUE ( S )  TO  BE  MODIFIED 

C  AND  MODIFY  THEM,  THEN  INVERT  THEM 

C 

CALL  VMQV < EV ( 1 ,  1),  1, EV< 1, 3),  1 ,  NH) 

JJ=11 

DO  444  JH=1,NH 

IF  ( IABS(KMAX( JH)-KREF).  GT.  3)  GO  TO  444 
EV< JH, 3 ) =EV ( NH,  3) 

PR ( JJi 6 ) =JH 
JJ=JJ+1 
444  CONTINUE 

CALL  VSDV(1.  ,  EV (1 ,  3 ) ,  1 , EV (1 , 4  > ,  1,  NH ) 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  CALCULATE  CONVENTIONAL  AND  MAXIMUM  LIKELIHOOD  BEAMFORMRING 

C 

DO  61  JB  =  1,NB£ 

LB=  <  JB-1 )*NHE+1 

CALL  WVMQV ( BSGW ( LB ) ,  1, BSQ< 1 >, 1, NHE) 

DO  62  JM=1,4 

CALL  VDOT ( P ( JB , JM ) , EV ( 1 , JM ) ,  1, BSG(1>,  1,  NH) 

62  CONTINUE 

DO  66  JH=1 ,  NH 

KMAX ( JH )  =  100.  *BSG< JH) 

IF  (JB.EQ.  NBE)  KMAX  (  JH )  =  100.  -#BSG  ( JH ) /QNRM 
66  CONTINUE 

KMAX (NHE )  =  100.  *BSG (NHE ) /QNRM 

IF  (JB.EQ. NBE)  KMAX <NHE)=100. *BSG(NHE) / (QNRM*QNRM> 

IF  ( ISSW(O) . GE.  0)  GO  TO  69 

WR ITE (6, 815)  JB, BGR ( JB  ) , ( KMAX ( JH >, JH=1 ,  NHE ),  JB 
815  FORMAT ( X,  13,  4X,  F10.  2,  4X,  2014,  4X,  16) 

69  CONTINUE 
61  CONTINUE 

CALL  VMIN (MIN, P ( 1 , 1), 1 , NB ) 

CALL  VSDV ( 1 .  ,P(1,2),1,P(1,2),1,NB) 

CALL  VSDV (  1 .  ,P(1,4),  1,P(1,4),  1,NB) 

QR=QS-P ( NBE,  2) 

C  WRITEC6, 845)  QS, P ( NBE, 2 ) , QR 

C  845  FORMAT (3F10.  3) 

IF  (ISSW(O).  LT.  0)  WRITEC6,  535) 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C  MAIN  PRINT  OUTPUT 

C 

DO  536  KA=1 , 2 

CALL  VSMYC100. , P ( 1 , KA ) , 1 , P ( 1 , KA ) , 1 , NB ) 

536  CONTINUE 
PR ( 1 , 6 >  =  IC 
PR ( 2,  6 ) =NA 
PRO,  6 )=JT 
PR  (4,  6 )  =  JF 
PR  <  5,  6 )  =F  ( JF ) 

PR  <  6,  6 )  =PQW 

PR ( 7 , 6 ) =GNRM* 100  /(I  *NH> 

PR  ( 8,  6  )  =GS*100. 

PR (9, 6)=GR*100. 

PR( 10, 6 ) =BGR ( KREF ) 

PR (20, 6 ) =MGTP 
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PR ( 1 9#  6 ) =C 

PR (18# 6)=F< JF>*10. 

PR ( 17 1  6)=A 
DO  876  JH=UJ, 16 
PR<  JH,  6 >=0. 

876  CONTINUE 

WR I TE( 6/ 538) 

538  FORMAT < 12X#  3(4X,  "  #  BEARING  CON  MLM  SPL  BE1  BE2  BQ"  >  ) 

DO  531  1*1. 3 
M=< 1-1 ) *2+1 
DO  533  0=1 i 20 
IJ=< 1-1 ) *20+ J 
DO  534  K=1  #  3 
NN  <  K )  =  <  K- 1 ) *60+ 1 J 

534  CONTINUE 

WRITE (6/ 9)  J, PR< J,  M),  PR< J,  M+l ) ,  <NN(K). BGR ( NN<  K) >, P(NN(K).  10), 
( P ( NN < K ) , L ) *  L= 1 , 2 ) »  <P(NN<K>, L ) >  L-6, 9), K*l.  3) 

9  FORMAT (314,  3  <  4X,  914) ) 

533  CONTINUE 
531  CONTINUE 

DO  537  KA»1,2 

CALL  VSMY ( .  01.  P<1.  KA).  1.  P<1.  KA).  1.  NB) 

537  CONTINUE 

DO  676  JB  =  1,NB 

PH=2.  5*ALN10/<1.  +  .  5*ABS  < P ( JB,  10)-BGR( JB ) ) ) 

P  ( UBi  9 )  *=EXP  <PH*P  <  JBi  6  )  /250.  ) 

P( JBi 10)=EXP (PH) 

676  CONTINUE 

WR ITE< 6 1  535) 

535  FORMAT ( 1H1 ) 

C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C  SELECT  REFERENCE  BEAM  DIRECTION  AND  DO  REFERENCE 
C  BEAM  NOISE  SUPPRESSION  . 

C 

C  NREF  -  REFERENCE  BEAM  NUMBER 

C 

NREF=KREF 

LREF=(NREF-1)*NHE+1 
DO  96  KA2*! i  2 

CALL  WVMGV ( BEW<  LREF# KA) #  1  #  BE < 1 # KA ) #  1,  NH) 

96  CONTINUE 

DO  97  JB=1 » NB 
LB=S(  JB-1  )*NHE+1 
DO  98  KA=1#  2 

CALL  WVMOV  <  BEW  <  LB  #  KA ) ,  1. B< 1. KA);  1.  NH) 

DO  99  KB=1 , 2 

CALL  VMPY(B<1,  KA ) ,  1.  BE(1.  KB).  1.  H(l.  KA.  KB).  1.  NH) 

99  CONTINUE 
98  CONTINUE 

CALL  VADD(H(1,  1,  1),  l.H(1.2.  2).  1.HC1.  1.  1);  1,  NH) 

CALL  VSUB  <  H(  1 ,  1;2),  l,H(i,2,  1);  1/H(1,2,  1),  1,  NH) 

DO  294  KM*=1;2 

UM=KM+4 

DO  293  KA=1 , 2 

CALL  VDOT (H(KA/  1. 2).  EV(1.  KM).  1.  HC1.  KA.  1 ).  1.  NH) 

293  CONTINUE 

CALL  VDOT  (P  ( UB;  JM),H<1.1,2).1.H<1.1.2>.1.2> 

294  CONTINUE 

97  CONTINUE 

SC  =  1.  /P (NREF#  1) 

CALL  VSMY (SC# P< 1.  5) ,  1. P< 1 , 5) ,  1.  NB ) 

CALL  VSMY (P (NREF#  2)#  P(l.  6),  1 #  P ( 1 , 6 > ,  1 ,  NB ) 

CALL  VMPY (P  < 1 #  6 ) ,  1 # P < 1 , 2 ) ,  1 , P ( 1 , 6 ) ,  1 ,  NB ) 

CALL  VMPY(P(1#  6),  1.  PCI.  2).  1.  P(l.  6).  1#NB) 

DO  292  KM=1,2 
JM=KM+4 

CALL  VSUB ( P ( 1 #  KM ) #  1 , P ( 1 , JM ) ,  1 . P ( 1 , UM > #  1 ,  NB > 

292  CONTINUE 
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ccccccccccccccccccccccccccccccccccccccccc 

c 

C  DO  REFERENCE  HYDROPHONE  NOISE  SUPPRESSION 

C 

CALL  VMOV ( 1 .  ,  0, EV  < 1 ,  1  > ,  1 ,  NH ) 

LREF=NB*NHE+i 
DO  601  KA»1,2 

CALL  WVMOV ( BEW ( LREF ,  KA),  1, BE(i, KA),  1,  NH) 

601  CONTINUE 

DO  602  JB®1,NB 
LB=< JB-1 )*NHE+1 
DO  603  KA=1 ,  2 

CALL  WVMOV  ( BEW  ( LB  ,  KA),  1,  B<  1,  KA),  1,  NH) 

DO  604  KB®  1,2 

CALL  VMPY  (B  ( 1 ,  KA),  1,  BE  (1,  KB),  1 ,  H  ( 1 ,  KA,  KB  ) ,  1,  NH) 

604  CONTINUE 
603  CONTINUE 

CALL  VADD (H C 1 ,  1,  1),  1,H(1,2,  2),  1,H(1,  1,  1),  1,  NH) 
CALL  VSUB  <  H( 1 ,  1 , 2 ) ,  1 , H ( 1 , 2,  1 ) ,  1 , H < 1 , 2,  1 ) ,  1 , NH ) 
DO  605  KM®!  ,  2 
JM=KM+6 
DO  606  KA=1 ,  2 

CALL  VDOT (H(KA,  1 , 2 > , EV ( 1 , KM) ,  1 , H< 1 , KA,  1 ) ,  1 ,  NH ) 

606  CONTINUE 

CALL  VDOT  (P  ( JB,  JM),H(1,1,2),1,H(1,1,2),1,2> 

605  CONTINUE 

602  CONTINUE 
SC»1.  /GS 

CALL  VSMY (SC, P < 1 , 7 ) ,  1 , P ( 1 , 7) ,  1 , NB ) 

SC— 1.  /QR 

CALL  VSMY (SC, P ( 1 , 8 ) , 1 , P < 1 , 8 ) , 1 , NB ) 

CALL  VSDV <  1 .  ,  P(  1, 2),  1,  P<  1,  2),  1 ,  NB  ) 

DO  607  KM=1 , 2 
JM=KM+6 

CALL  VSUB  <  P  ( 1 ,  KM),  1 ,  P  <  1 ,  JM  > ,  1»  P<1«  JM>«  1,  NB) 

607  CONTINUE 

CALL  VSDV ( 1 . , P ( 1 , 2 ) , 1 , P ( 1 , 2 ) , 1 , NB ) 

CALL  VSDV (  1 .  ,  P<1, 8),  1,  P(l,  8),  1 ,  NB ) 

IF  < ISSWC 1).  GE.  0)  GO  TO  609 
DO  281  JB  =  1,NB 

WRITE <6, 816)  JB, BGR ( JB ) ,  <P< JB, JM), JM=1, 10), JB 
816  FORMAT ( X,  13,  4X,  F10.  2,  4X,  10F10.  4,  4X,  16) 

281  CONTINUE 

WRITE <6,  535) 

609  CONTINUE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C  MOVE  P< JB, JM), JB®1,  NB  AND  JM=1,2  TO  DISC 

C 

JT1=JT 

IF(JT.  GT.  16)  JT1®JT-16 
LD=< JF-1 >*2 
ISEC=< JT1-1 >*6 
ITRK1=LD*NP 
DO  400  JD®1,NP 
ITRK=ITRKl+< JD-1 )*2 
IF(JT. GT. 16)  ITRK=ITRK+1 
C  DO  500  ID=1 , NB 

C  RB  ( ID )  ®P  ( ID,  JD ) 

C  500  CONTINUE 

C  WRITE<6,  599)  <P ( JB, JD) ,  JB=i,  NBE) 

CALL  INTER (P ( 1 ,  JD) ,  RB < 1 ) , F< JF) , NB ) 

CALL  VMAX  <NMAX, RB ( 1 ) ,  1 , NB ) 

RB ( 182 ) ®RB  <  NMAX ) 

RB ( 183 ) ®POW 
C  NBC=NB+2 

C  WRITE(6, 599)  <RB ( JB ) , JB=1 , NBC ) 

C  599  FORMAT < 10F8.  3) 

CALL  EXEC ( 2, 14, IB, 384, ITRK, ISEC) 

400  CONTINUE 
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0701  c 

0702  21  CONTINUE 

0703  C 

0704  C  END  OF  FREQUENCY  LOOP 

0705  C 

0706  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
0707  C 
0708  C 
0709  C 

0710  C  END  OF  MAIN  LOOP 

0711  C 

0712  191  CONTINUE 

0713  C 

0714  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
0715  C 

0716  C  ALL  DATA  ARE  PROCESSED  AND  THE  RESULTS  ARE  ON  DISC 

0717  C 

0718  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
0719  C 

0720  C  END  OF  MAIN  PROGRAM 

0721  C 

0722  IMERD»LREC(1) 

0723  WRITE (1*9191)  IMERD 

0724  9191  FORMAT ( /  *  "  LAST  BLOCK  :  ",  16) 

0725  C 
0726  C 

0727  END 

FTN4  COMPILER:  HP92060-16092  REV.  2001  (791101) 

**  NO  WARNINGS  **  NO  ERRORS  **  PROGRAM  *  10953  COMMON  -  00000 
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E. 2  SUBROUTINE  EIGEN 

The  EIGEN  subroutine  calculates  the  eigenvalues  and  eigenvectors  of  a 
Hermitian  matrix.  A  listing  of  the  subroutine  follows. 


0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 

0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 

0056 

0057 

0058 

0059 

0060 

0061 

0062 

0063 

0064 

0065 

0066 


FTN4,  L 

SUBROUTINE  EIGEN( ARE# AIM#  ERE#  EIM#  ISIZE#  EIGV, N> 
DIMENSION  ARE ( ISIZE)#  AIM (ISIZE ), ERE ( IS I ZE ) #  E I M ( ISI ZE ) 
DIMENSION  EIGV ( 20 ) #  H ( 20#  4 ) 

INTEGER  P,Q#NN(20) 

C 

C  DATA  FOR  BREAKING  OFF 

C 

NNN=N 

C 

IOW=N 

C 

EPS=.  IE— 9 
C 

C  D  UNITY  MATRIX  ERE# EIM# FOR  EIGENVECTORS 

C 

DO  1  1  =  1#  N 
DO  1  U=1#N 
IJ»(I-1)*N+J 
JI=(J-1 >*N+I 
ERE< I  J)=0. 

E IM ( I J )  =0. 

IF  (I.GT.  J)  ARE(IJ)=ARE<JI)*ARE(JI>+AIM<JI)*AIM<JI> 

1  CONTINUE 
DO  2  1=1# N 
K9=( 1-1 )*N+I 
ERE(K9)=1. 

2  CONTINUE 
C 

C  BEGIN  OF  ITERATION 

C 

DO  100  IT=1 # 1000 
C 

C  COMPUTATION  OF  SQUARE  SUM  OF  NONDIAGONAL  ELEMENTS 

C  AND  DETERMINATION  OF  MAX  NONDIAGONAL  ELEMENT 

C 

C  SW=0. 

H(  1 1  1 )  =0, 

DO  101  1=2# N 
NEL=I-1 
IU=( 1-1 >*N+1 

CALL  UMAX ( NN ( I  ) #  ARE ( I J ) #  1 #  NEL  > 

H ( I # 1 )=ARE< IJ-1+NN< I > > 

101  CONTINUE 

CALL  VMAX  ( Q#  H(l#  1  )#  1#  N) 

P=NN(G) 

SW=SQRT(H(Q,  1 ) ) 

NNN= I T 

IF(SW. LT.  EPS)  GOTO  899 
C 

C  PIVOT-ELEMENT  FOUND  COMPUTATION  ROTATION  PARAMETERS 

C 

NP=N*P 

NG=N*G 

NPP=NP-N+P 

NPG=NP-N+Q 

NGP=NG-N+P 

NGG=NG-N+G 

APP=ARE ( NPP ) 

AGG=ARE(NGG) 

APG=ARE<NPG) 

BPG=AIM (NPG) 

RO=SQRT ( APG*APG+BPQ*BPG ) 

R0D=2.  *R0 

C2X» ( AGG-APP ) /ROD 

VZ=-1. 
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0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 

0081 

0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 

OlOO 

0101 

0102 

0103 

0104 

0105 

0106 

0107 

0108 

0109 

0110 

0111 

01 12 

0113 

0114 

0115 

0116 

0117 

0118 

0119 

0120 

0121 

0122 

0123 

0124 

0125 

0126 

0127 

0128 

0129 

0130 

0131 

0132 

0133 

0134 

0135 

0136 

0137 

0138 


SACLANTCEN  SM-158 


IF < C2X.  GE.  0 )  VZ*1 
CX=C2X+VZ*SGRT(C2X*C2X+1  ) 

C2X=CX*CX 
CC=C2X/<1.  +C2X) 

C=SQRT (CC ) 

S=C/CX 

SS=S*S 

SOS*C 

SR=APG*S/RO 

SI=~BPG*S/RO 

ROTATION 

NEL=P-1 

IP*P 

IG=Q 

CALL  VSMY ( SI , AIM< IG ) ,  N,  H (1 ,  1 >,  1, NEL) 

CALL  VSMY ( —SR, AIM(IG), N, HC1, 2), 1, NEL) 

CALL  VSMY (SI , AIM(IP)i N, H(l, 3). 1, NEL) 

CALL  VSMY ( SR i  AIM(IP),N,  H(l,4>,  1,  NEL) 

CALL  VPIVC-SR,  ARE(IG),  N,  HC1,  l),liH(Ll)t  1  ,  NEL ) 
CALL  VPIVC-SI,  ARE<  IG) ,  Ni  H(l,  2) ,  1,H(1,2>,  1,  NEL) 
CALL  VPIV(SRi  ARE (  IP),N,  H<1,3),1,HU,3),1,  NEL) 
CALL  VP IV ( -SI , ARE (IP), N, H ( 1 , 4 ) ,  1 , H ( 1 , 4 ) ,  1 , NEL) 
CALL  VPIV(C, ARE< IP ) , N, H ( 1 ,  1 ) ,  1 ,  ARE( IP), N, NEL) 
CALL  VPIVCC, AIM(IP), N, H<1, 2),  1,  AIM(IP),  N,  NEL) 
CALL  VPIV<C, ARE( IG), N, H< 1, 3),  1,  ARE( IG),  N,  NEL) 
CALL  VPIVCCi  AIM(  IG),  N,  H<  1,  4),  1,  AIM<  IG),  N,  NEL) 
JP-NP-N+ 1 
JG=NQ-N+1 

CALL  VMPY< ARE( IP),  N,  ARE( IP),  N,  AREC JP),  1, NEL) 
CALL  VMPY< AIM< IP ) ,  N,  AIM(IP), N, AIM(JP),  1, NEL) 
CALL  VMPY(ARE( IG),  N,  ARE( IG),  N,  ARE( JG),  1, NEL) 
CALL  VMPY ( A IM < IG ) , N, A IM ( IG ) , N,  AIM ( JG ) ,  1, NEL) 
CALL  VADD ( ARE (JP),  1, AIM( JP),  1, ARE( JP),  1, NEL) 
CALL  VADD(ARE( JG), 1, AIM( JG), 1, ARE( JG), 1, NEL) 
NEL=Q-P-1 
IP-NPP+1 
IG-NP+G 

CALL  VSMY (SI,  AIM (IQ),  N,  H(l,  1),  1, NEL) 

CALL  VSMY  (SR,  AIM  ( IQ ) ,  N,  HU,  2),  1,  NEL) 

CALL  VSMY  ( -SI,  AIM  ( IP),  1,  H(l,  3),  1,  NEL) 

CALL  VSMY  ( -SR ,  A I M  (  IP),  1,H(1,4>,  1,  NEL) 

CALL  VPIVC-SR,  ARE(  IG),  N,  HU,  1),1,HU,1),  1,  NEL) 
CALL  VPIV(SI,  ARE ( IG ) ,  N,  H (1 ,  2 ) ,  1 , H (1 , 2 > ,  1 , NEL) 
CALL  VP  IV  ( SR,  ARE  ( IP  ) ,  1,H<1,3>,  1,HU,3>,  1,  NEL) 
CALL  VPIV(-SI,  ARE ( IP ) ,  1, H(l, 4),  1, H(l, 4),  1,  NEL) 
CALL  VPIVCC, ARE ( IP), 1, H(l, 1 ), 1, ARE (IP), 1, NEL) 
CALL  VPIV(C,  AIM(  IP),  1,H(1,2>,  1,  AIM  (IP),  1 ,  NEL  ) 
CALL  VPIV(C,  ARE  ( IG ) ,  N,  HU,  3),  1,  ARE  (IQ),  N,  NEL) 
CALL  VPIV(C,  AIM  (  IG ) ,  N,  H(l,  4),  1,  AIM  (IQ),  N,  NEL) 
JP=NP+P 
JG=NGP+1 

CALL  VMPY (ARE ( IP ) ,  1, ARE(IP),  1, ARE(JP), N,  NEL) 
CALL  VMPY (AIM (IP),  1, AIM (IP),  1, AIM (JP), N, NEL) 
CALL  VMPY ( ARE ( IQ),  N,  AREC IG),  N,  ARE( JG),  1, NEL) 
CALL  VMPY (AIM(IQ),  N,  AIM( IG) ,  N,  AIM( JG ) ,  1, NEL) 
CALL  VADD ( ARE (JP),  N,  AIM( JP),  N,  ARE(JP),  N,  NEL) 
CALL  VADD(ARE( JG),  1, AIM(JG),  1, ARE( JG),  1, NEL) 
NEL=N-G 
IP=NPG+1 
IG=NGG+1 

CALL  VSMY ( -SI,  AIM(  IG) ,  1,H(1,  1),  1 ,  NEL ) 

CALL  VSMY (—SR,  AIM(IG),  1, H(l, 2),  1, NEL) 

CALL  VSMY (-SI ,  AIM(IP),  1, H(l, 3),  1, NEL) 

CALL  VSMY  (SR,  AIM(  IP),  1,H(1,4),  1 ,  NEL) 

CALL  VP IV(-SR,  ARE(IQ),  1 , H (1 ,  1 ) ,  1 , H(1 ,  1 ) ,  1 , NEL) 
CALL  VPIV(SI,  ARE ( IG > ,  1,H(1,2),1,HU,2),1,  NEL) 
CALL  VPIV(SR,  ARE  ( IP  ) ,  l,HU,3),l,H(l,3),i,  NEL) 
CALL  VPIVCSI,  ARE(IP),  l,H(i,4),l,HU,4),l,  NEL) 
CALL  VPIVCC,  ARE(IP),  1,H(1,  1),  1,  AREC  IP),  1 ,  NEL ) 
CALL  VPIVCC,  AIM(  IP),  i,HCl,2>,  1,  AIM  (IP),  1 ,  NEL ) 
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0139 

0140 

0141 

0142 

0143 

0144 

0145 

0146 

0147 

0148 

0149 

0150 

0151 

0152 

0153 

0154 

0155 

0156 

0157 

0158 

0159 

0160 

0161 

0162 

0163 

0164 

0165 

0166 

0167 

0168 

0169 

0170 

0171 

0172 

0173 

0174 

0175 

0176 

0177 

0178 

0179 

0180 

0181 

0182 

0183 

0184 

0185 

0186 

0187 

0188 

0189 

0190 

0191 

0192 

0193 

0194 

0195 

0196 

0197 

0198 

0199 

0200 

0201 

0202 

0203 

0204 

0205 

0206 

0207 

0208 

0209 


CALL  VPIVCC,  ARE(IG),  1,H<1,3),  1,  ARE <  IQ),  1,NEL> 

CALL  VPIV(C,  AIM  (IQ),  1,H(1,4),  1,  AIM(IG),  1#  NEL ) 

JP=NG+P 

JG=NG+G 

CALL  VMPY(ARE(IP>,  1, ARE(IP),  1, ARE(UP),  N, NEL) 

CALL  VMPY ( AIM ( IP ) ,  1, AIM (IP),  1#  AIM(UP  >, N, NEL) 

CALL  VMPY(ARE(IG>,  1, ARE(IQ),  1,  ARE(JG), N, NEL) 

CALL  VMPY(AIMUQ),  1,  AIM(IG),  1,  AIM(JG),  N,  NEL) 

CALL  VADD ( ARE ( JP  ) ,  N,  AIM(JP),  N,  ARE ( JP  ) ,  N, NEL) 

CALL  VADD ( ARE ( JG ) , N,  AIM(JG),  N,  ARE( JG),  N,  NEL) 

C 

C  E.  B.  L.  18/02/82 
C 

ARE(NPP)=APP*CC+AGG*SS-ROD*SC 

ARE(NGQ>=APP*SS+AGG*CC+ROD*SC 

ARE(NPG>=0. 

AIM(NPG)=0. 

ARE(NGP)=0. 

C 

C  COMPUTATION  OF  EIGENVECTORS 

C 

IP=NP-N+1 

IG=NG-N+1 

CALL  VSMYCSI, EIM(IQ),  l,H(Ll)iLN) 

CALL  VSMY <  —SR,  EIM<  IG),1,H(1,2),1,N> 

CALL  VSMY  (SI,  EIM<  IP),1,H(1,3),1,N) 

CALL  VSMY  (SR,  EIM(  IP),  1,H(1,4),1,N) 

CALL  VPIVC-SR,  ERE(  IG),1,H<1,1),1,H<1,1>,1,N) 

CALL  VPIVC-SI, ERE(IG),  1,  H( 1,  2) ,  1,  H( 1, 2) ,  1,  N) 

CALL  VP IV < SR i  ERE ( IP ) ,  1,  H<  1,  3),  1,  H<  1,  3) ,  1,  N) 

CALL  VPIV(-SI,  ERE( IP ) ,  1 ,  H( 1#  4) ,  1 ,  H< 1  *  4 ) ,  1 ,  N) 

CALL  VPIV(C, ERE(IP),  1, H(l,  1).  1,  ERE(IP),  1, N) 

CALL  VPIV(C,  EIM(IP),  1,  H(l,  2),  l.EIM(IP),  1,  N) 

CALL  VPIV(C,  ERE<  IG )  ,  i,H(l,3),  1 ,  ERE(  IG ) ,  1,  N) 

CALL  VPIV(C, EIM( IG ) ,  1 ,  H< 1 , 4 ) ,  1 ,  EIMC IG),  1, N) 

1006  FORMAT (2X,  'ARE=',  4<F10.  3,  3X) .  5X,  'AIM=',  4(F10.  3,  3X)  ) 

1007  FORMAT  <2X,  'ERE-  ' ,  4<F10.  3,  3X),  5X,  /EIM=/,  4(F10.  3,  3X )  ) 

C 

C  ROTATION  END 

100  CONTINUE 
C 

C  200  ITERATIONS  PERFORMED  SET  N=-N 

C 

NNN=-N 

C  WRITEC1,  1002)  NNN 

C 1002  FORMAT  <  10X,  'NNN=-N=  ' ,  13  > 

C 

C  SORT  OF  EIGENVALUES  WITH  INCREASING  ABS  VALUES 

C 

899  LL=1 

900  AK=-1.E04 

DO  901  1  =  1,  N 
K2=(I-1)*N+I 

901  AK=AMAX1(AK,  ARE(K2  > > 

DO  902  1  =  1,  N 

Kl=< 1-1 )*N+I 

IF(AK.  EG. ARE(Kl))  GOTO  903 

902  CONTINUE 
C 

C  FOUND  STORE  IN  EIGV  MARC  POSITION 

C 

903  NN(LL)=I 
EIGV(LL)=AK 
LL=LL+1 
Kll=< I— 1 )*N+I 
ARE(K1 1 )=-l.  EOS 
IF(LL-N)  900,900,904 

C 

C  EIGENVALUES  ARE  SORTED  IN  EIGV.  SORT  OF  COLUMNS 

C  OF  ERE,  EIM  IN  THE  SAME  ORDER  AND  STORE  IN  ARE, AIM. 

C 
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0210 

0211 

0212 

0213 

0214 

0215 

0216 

0217 

0218 

0219 

0220 

0221 

0222 

0223 

0224 

0225 

0226 

0227 

0228 

0229 


904  DO  905  1  =  1,  N 
K0=< 1-1 )*N+1 
K10=(NN!I)-1)*N+1 

CALL  VMOV (ERE!K10 ) ,  1,  ARE(KO),  1,  N) 

CALL  VSMY(~1.  , EIM(KIO) ,  1, AIM(KO),  1,  N) 

905  CONTINUE 
N=NNN 
ICL=1 5414B 
IBL=3407B 

CALL  EXEC <2, 1001B, ICL, 1) 

CALL  EXEC<2» 1001B, IBL* 1) 

WRITE! 1, 1003)  N 

1003  FORMAT! 1 OX, 'OUT  OF  EIGEN  N=NNN= 13) 
WRITE! 1 , 2004 )  (EIGVCI), 1=1, IOW) 

2004  FORMAT  C  5X,  'EIGV=  ',F15.  5) 


RETURN  BREAKING  OFF 


RETURN 
END 

FTN4  COMPILER:  HP92060-16092  REV. 
**  NO  WARNINGS  **  NO  ERRORS  ** 


2001  (791101) 
PROGRAM  =  02935 


COMMON  =  00000 
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E. 3  ASSEMBLER  PROGRAM  FOR  DATA  TRANSFER 

The  following  is  a  listing  of  a  program  for  transferring  data  from  magnetic 
tape  into  a  part  of  the  off-line  HP  computer  memory  called  "EMA11. 


****  COM, V.  DUARTE  REV.  00  05/03/81  **** 


COM 


V.  DUARTE 


MARCH  81 


REV  00 


0001  ASMB, B#  L, T 

0002  00000  NAM  LEMA1 #  7 

0003  ENT  LEMA1 

0004  EXT  .  EMAP#  .  ENTR,  EXEC 

0005  EXT  RBE 

0006* 

0007****************************************************************** 

0008* 

0009* 

0010* 

0011* 

0012* 

0013* 

0014* 

0015* 

0016* 

0017* 

0018* 

0019* 

0020* 

0021* 

0022* 

0023* 

0024* 

0025* 

0026****************************************************************** 


CALLING  SEQUENCE 


CALL  LEMA1 (LU, I ERR ) 

LU  -  MAGNETIC  TARE  LOGIC  UNIT 
I ERR  -  RETURN  IF  ERROR 

PURPOSE  : 

ROUTINE  TO  TRANSFER  DATA  FROM  MAGNETIC  TAPE  INTO 
AN  “EMA"  ARRAY. “EMA"  ARRAY  SHOULD  BE  DECLARED  AS  THE 

FIRST  OF  THE  “EMA”  AREA. 


0027* 

0028 

00000 

000000 

LU 

NOP 

MT  LOGIC  UNIT 

0029 

00001 

000000 

I  ERR 

NOP 

ERROR  RETURN 

0030 

00002 

000000 

LEMA1 

NOP 

ROUTINE  ENTRY  POINT 

0031 

00003 

01 6002 X 

USB 

.  ENTR 

GET  INPUT  PARAMETER 

0032 

00004 

OOOOOOR 

DEF 

LU 

0033 

00005 

0 1 600 1 X 

JSB 

.  EMAP 

GET  USER  MAP  ADDRESS 

0034 

00006 

000012R 

DEF 

*4-4 

0035 

00007 

000004X 

DEF 

RBE 

OF  CORRESPONDENT  EMA 

0036 

00010 

000032R 

DEF 

TABLE 

AREA  ARRAY 

0037 

00011 

000037R 

DEF 

INDEX 

0038 

00012 

026030R 

RTN 

JMP 

ERROR 

SEND  ERROR  MESSAGE  TO  MAIN 

0039 

00013 

07604 1R 

STB 

BUF 

SAVE  ADDRESS 

0040 

00014 

162000R 

LDA 

LU,  I 

PREPARE  CONTROL  WORD 

0041 

00015 

042043R 

ADA 

=*B  100 

FOR  EXEC  CALL  AS  BINARY  READ 

0042 

00016 

072000R 

STA 

LU 

0043 

00017 

016003X 

JSB 

EXEC 

READ  DATA  FROM  MT 

0044 

00020 

00002 5R 

DEF 

*•+•5 

0045 

00021 

000040R 

DEF 

CODE 

0046 

00022 

OOOOOOR 

DEF 

LU 

0047 

00023 

100041R 

DEF 

BUF#  I 

0048 

00024 

000042R 

DEF 

NWORD 

0049 

00025 

002400 

CLA 

0050 

00026 

172001R 

END  I 

STA 

IERR, I 

0051 

00027 

126002R 

JMP 

LEMA1 # I 

0052 

00030 

003400 

ERROR 

CCA 

0053 

00031 

026026R 

JMP 

END  I 

0054 

00032 

000001 

TABLE 

DEC 

1 

0055 

00033 

177777 

DEC 

-1 

0056 

00034 

000001 

DEC 

1 

0057 

00035 

000000 

DEC 

0 

0058 

00036 

OOOOOO 

DEC 

0 

0059 

00037 

000001 

INDEX 

DEC 

1 

0060 

00040 

000001 

CODE 

DEC 

1 

0061 

00041 

OOOOOO 

BUF 

NOP 

0062 

00042 

001457 

NWORD 

DEC 

815 

00043 

000100 

0063 

END 

**  NO  ERRORS  *TOTAL  **RTE  ASMB  92067- 

16011** 
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E. 4  PROGRAM  FOR  ON-LINE  PROCESSING 

The  purpose  of  the  on-line  processing  is  described  in  Chapter  2.  The 
following  is  a  listing  of  the  program. 


&SPM8M  T-00004  IS  ON  CR00013  USING  00036  BLK5  R“0000 


0001 

0002 

0003 

0004 

0003 

0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0013 

0016 

0017 

0018 

0019 

0020 

0021 

0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 

0056 

0057 

0058 

0059 

0060 

0061 

0062 

0063 


FTN4# L 

*EMA< IBEMA, 5) 

PROGRAM  SUPVR 

DIMENSION  SF < 20 ) »  NDXSF ( 20 ) 

DIMENSION  I HEAD < 15), IB1 <4) , RFG<22> 

C  DIMENSION  ITIME(IO) 

DIMENSION  PQWR ( 184  > , NAM1 (3) ,  IFQ<2> 

DIMENSION  IDCB(IO) 

COMMON  /IBEMA/  IBB<204S0) 

COMPLEX  CX 

COMMON  LUMT, IHEAD, CX ( 400 ) 

EQUIVALENCE  ( IB1 ( 1 ) , CX < 1 ) ) »  (RFG< 1 ) #  IB1<4)> 

EQUIVALENCE  ( IFQ ( 1 ) , SF< 1 ) ) 

EQUIVALENCE  ( AFRST,  IDCB(3>  >,  (ALAST,  IDCB(5) ),  (DBM IN,  IDCB(7) /, 

*  (DBMAX, IDCB (9) ) *  (NPTS, POWR <1 ) ), (PMAX, P0WR<2) ) 

DATA  IDCB  <  1  )  i  DBM  IN,  DBMAX /O,  -30.  0,0.  0/ 

DATA  NAM1/2HBF, 2HPL,  2HT  / 

DATA  IBUS2, IBUS3/12B, 192/ 

DATA  IEX,  IST/2HEXi  2HST/ 

DO  15  K*l, 15 
15  IHEAD  (K)  ==0 

LUMT»10 

10  WRITE  < 1 1  100 ) 

100  FORMAT ( 'NUMBER  OF  SAMPLES ( POWER  OF  2)?') 

READC 1,  * )  NS 
1*1 

20  IF  ( I -NS )  25,40,30 

25  1=1+1 

GOTO  20 

30  WRITE<  1,200) 

200  FORMAT < 'NUMBER  OF  SAMPLES  MUST  BE  A  POWER  OF  2'/ 

*  'PLEASE  ENTER  A  NEW  VALUE') 

GOTO  10 

40  RNS=NS 

WRITE< 1,300) 

300  FORMAT ( 'NUMBER  OF  HYDROPHONES,  SP AC ING ( METRES )?' ) 

READd,*)  NH,  D 

IF<NH  .  LE.  0)  GOTO  10 

RNH«NH 

BS»RNS*RNH 

IF ( BS  .LE.  40960.0)  GOTO  45 
WRITEd,  400) 

400  FORMAT < 'CAPACITY  OF  MAP  MEMORY  BUS  2  EXCEEDED'/ 

*  'REDUCE  NUMBER  OF  SAMPLES  AND/OR  HYDROPHONES') 

GOTO  40 

45  WRITEd,  450) 

450  FORMAT < 'NUMBER  OF  BEAMS,  INITIAL  ANGLE  «c  FINAL  ANGLE < DEGREES )?' ) 

READd,*)  NBMS,  AFRST,  ALAST 
C 

C  COMPUTE  SPACE  REQUIRED  FOR  CO-VARIANCE  MATRICES 

C  AND  STEERING  VECTORS  FOR  A  SINGLE  FREQUENCY 

C 

SP=2. 0*RNH*(NH+NBMS)-2*NBMS 
IF (SP+BS  .  LE.  40960.0)  GOTO  46 
WRITEd,  470) 

470  FORMAT < 'NO  SPACE  FOR  CO-VARIANCE  MATRICES  &  STEERING  VECTORS ' ) 
CALL  EXEC <6 ) 

46  WRITEd,  460) 

460  FORMAT< 'SAMPLING  FREQUENCY?') 

READd,*)  SFQ 
DF=SFG/NS 
FQMAX=DF*< NS/2-1) 

50  M~40960,  O-BS 
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0064 

0065 

0066 

0067 

0068 

0069 

0070 

007  i 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 

0081 

0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 

0100 

0101 

0102 

0103 

0104 

0105 

0106 

0107 

0108 

0109 

0110 

0111 

0112 

0113 

0114 

0115 

0116 

0117 

0118 

0119 

0120 

0121 

0122 

0123 

0124 

0125 

0126 

0127 

0128 

0129 

0130 

0131 

0132 

0133 

0134 

0135 
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NN“2*NH*NH 

C  NN=SIZE  OF  EACH  CO- VARIANCE  MATRIX 

L=20 

60  WR ITE< 1 1  600 ) 

600  FORMAT < 'HOW  MANY  FREQUENCIES  DO  YOU  WISH  TO  SELECT? ' ) 

READ<1,*)  NF 
IF<NF  .  LE.  0)  GOTO  40 
IF(NF  .LE.  L)  GOTO  70 
WRITEU,  700) 

700  FORMAT( 'MAXIMUM  NUMBER  OF  FRQUENCIES  EXCEEDED'/ 

*  'PLEASE  ENTER  A  NEW  VALUE') 

GOTO  60 

70  RNF“NF 

WRITE ( 1/ 800 )  DF ,  FQMAX 

800  FORMAT ( 'FREQUENCY  RESOLUTION  =  ',FB.  2,  '  HZ'/ 

*  'MAXIMUM  FREQUENCY  IS  ',F6.  1,  '  HZ'/ 

*  'INPUT  SELECTED  FREQUENCIES  USING  FREE  FIELD  FORMAT'/ 

*  'USE  SPACE  OR  COMMA  TO  SEPARATE  DATA  ITEMS'/ 

*  'TERMINATE  EACH  LINE < EXCEPT  THE  LAST  ONE)  BY  /'/ 

*  'FIRST  FREQUENCY  SPECIFIED  SHOULD  BE  THE  ONE'/ 

*  'FOR  IMMEDIATE  PROCESSING  AND  DISPLAY') 

READ ( 1/ * )  <  SF  < I ) »  1  =  1, NF) 

DO  72  1*1, NF 

I F < SF < I )  .  GT.  FQMAX)  GOTO  73 
K=SF  < I ) /DF+O.  5 
SF  < I ) =K*DF 
NDXSF(I)=K+1 

72  CONTINUE 
GOTO  74 

73  WRITEU,  480) 

480  FORMAT < 'MAXIMUM  FREQUENCY  EXCEEDED'/ 

*  'PLEASE  SELECT  A  NEW  SET') 

GOTO  70 

74  WRITE< 1,500) 

500  FORMAT< 'SPECIFY  PROCESSING  TYPE '/ ' 1=CBF '/ '2=CBS ' / '3=MLM ' / '4=MLS ' > 

READ  < 1,  *  >  I TYPE 

IF< ITYPE  .  LT.  1  .OR.  I TYPE  .  GT.  4)  GOTO  74 

WRITEU,  1000) 

1000  FORMAT < 'NUMBER  OF  DATA  ACQUISITIONS?') 

READ  U ,  * )  NA 
C 

C  PASS  FIRST  BUFFER  TO  PLOT  PROGRAM 

C  THEN  SCHEDULE  PLOT  PROGRAM 

C 

NPTS=NBMS 
I LEN=2* ( NBMS+3 ) 

ICLS=0 

CALL  EXEC (20, 0, IDCB, 10, 0, 0, ICLS) 

CALL  EXECUO,  NAM1 ,  ICLS,  ITYPE,  IFQ(l),  IFG<2),  ILEN ) 

C 

WRITE ( 1 # 75) 

75  FORMAT ( 'RUN-ID< 1-8  CHARS)?') 

READ( 1 , 76 )  < IHEAD (K ) ,  K=10,  13) 

76  FORMAT <4A2) 

D  WRITE<6,  950) 

D950  FORMAT ( 1H  , 'FREQUENCY  INDEX') 

D  DO  75  1*1 , NF 

D  WRITEC6, 980)  SF< I ) , NDXSF < I ) 

D980  FORMAT < 1H  ,  F9  1,  19) 

D75  CONTINUE 

WRITEU,  1100) 

1100  FORMAT ('TO  TERMINATE  THE  PROGRAM  HIT  ANY  KEY'/ 

*  'THEN  GIVE  THE  COMMAND  BR,SUPVR') 

C 

C  THE  DATA  ACQUISITION  PROGRAM  IS  EXPECTED  TO  PLACE  INFORMATION 

C  INTO  THE  FIRST  NS*NH  W0RDS<32  BITS  EACH)  OF  MEMORY  BUS  2 

C  SPACE  FOR  COVARIANCE  MATRIX  IS  ALLOCATED  AT  HIGH  END  OF  BUS  2 

C  BUT  LAST  WORD  IS  LEFT  FREE  TO  AVOID  BUFFER  ADDRESS  OUTSIDE  LIMITS 

C 

BSCV=NN 

BACV=81918.  0-2.  0*BSCV 
CALL  MPOPN (6 ) 
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0136 

0137 

0138 

0139 

0140 

0141 

0142 

0143 

0144 

0145 

0146 

0147 

0148 

0149 

0150 

0151 

0152 

0153 

0154 

0155 

0156 

0157 

0158 

0159 

0160 

0161 

0162 

0163 

0164 

0165 

0166 

0167 

0168 

0169 

0170 

0171 

0172 

0173 

0174 

0175 

0176 

0177 

0178 

0179 

0180 

0181 

0182 

0183 

0184 

0185 

0186 

0187 

0188 

0189 

0190 

0191 

0192 

0193 

0194 

0195 

0196 

0197 

0198 

0199 

0200 

0201 

0202 

0203 

0204 

0205 

0206 

0207 

0208 

0209 
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c 

C  SET  UP  THE  FOLLOWING  CONSTANTS  IN  THE  SCALAR  TABLE 

C  50  =  RNH*RNH, USED  IN  MLM  &  MLS  CALCULATIONS 

C  51  *  1. 0/ (RNH-1.  0) , USED  IN  CALCULATING  AVERAGE  POWER 

C  52-53  *  CMPLXU.  0,  0.  0),  USED  BY  CSDOT 

C  54  =  1.  0/(NA*NS*NS), USED  TO  NORMALISE  CO-VARIANCE  MATRIX 

C  55  =  100. 0/ USED  TO  SCALE  RESULTS  BEFORE  PRINTING 

C 

CX(  1 )®CMPLX(RNH*RNH,  1.  0/(RNH-l.  0) > 

CX  <  2 )  =1.  O 

CX  <3 ) “CMPLX ( 1 .  0/(NA*RNS*RNS),  100.  0) 

CALL  MPWST< 50i CX* 6»  1) 

C 

C  CALCULATION  OF  STEERING  VECTORS 

C 

PI=*3. 14159265 

A I NC  ®  < ALAST- AFRST ) / (NB MS- 1 ) 

A*AFRST 

C 

BSSV®2. 0* (RNH-l .  0)*NBMS 
BASV=BACV-2.  0*BSSV 
D  WRITE (6, 86)  BASV 

D86  FORMAT < '  STEERING  VECTORS  START  AT  BUS2  ADDRESS:  ',  F6.  1) 
LROW»NBMS 
ROWL»LROW 
C 

C  CONFIGURE  FIRST  ROW  AS  A  COMPLEX  VECTOR 

C  AND  FILL  WITH  THE  CONSTANT  (1.0,0.  0) 

C 

CALL  MPCLB  <21 # BASV,  ROWL,  1, NH-1, 0) 

DO  210  1*1 i LROW 
CX(X>»1.0 
210  CONTINUE 

CALL  MPWDB (21 , CX, !• 1) 

C 

C  CONFIGURE  SECOND  ROW  AND  INITIALIZE  TO 

C  EXP ( J*2*P I *F*D/C*CQS ( ANGLE ) > 

C 

BA«BASV+4. 0 

CALL  MPCLB <21, BA, ROWL,  1,  NH-1, 0) 

K=1 

C0NST»2.  0*PI*D*SF(1>/1500. 0 
DO  220  J»1,NBMS 
F=C0NST*C0S(PI*A/180.  0) 

CX  (K ) “CMPLX ( COS ( F ) i SIN(F> ) 

K=K+1 
A^A+AINC 
220  CONTINUE 

CALL  MPWDB ( 21 , CX,  1,  1 ) 

C 

C  SET  UP  REMAINING  ROWS  OF  MATRIX 

C  R0W(K)=R0W(2)*R0W(K-1 ) 

C 

CALL  MPEQB<22,  21) 

LST=NH-1 
DO  240  I<3,  LST 
BA=xBA+4.  0 

CALL  MPCLB (23, BA,  ROWL,  1, NH-1, 0) 

CALL  CVMUL(23,  0, 21, 0,  22) 

CALL  MPEQB(22,  23) 

240  CONTINUE 
C 

C  SET  UP  COSINE  TABLE  ON  BUS  3  FOR  USE  BY  FFT  FUNCTION 

C 

CTS=NS/2 

CALL  MPCLB(30,  0.0,  CTS,  0,  1,  0) 

CALL  VCOST ( 30,  CTS ) 

C 

C  SET  LST  TO  NO.  OF  HYDROPHONES- 1 

C  SET  RNK«NO.  OF  16-BIT  HALFWORDS  PER  COMPLETE  COVARIANCE  MATRIX 

C  SET  SNK»NO.  OF  16-BIT  HALFWORDS  PER  EXTRACTED  COVARIANCE  MATRIX 

C  INITIALIZE  BAC  TO  START  ADDRESS  OF  FIRST  COVARIANCE  MATRIX 

C  INITIALIZE  DAC  TO  START  ADDRESS  OF  FIRST  EXTRACTED  MATRIX 

C  INITIALIZE  BAG  TO  START  ADDRESS  OF  FIRST  Q-VECTOR ( COLUMN) 
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0210 

0211 

0212 

0213 

0214 

0215 

0216 

0217 

0218 

0219 

0220 

0221 

0222 

0223 

0224 

0225 

0226 

0227 

0228 

0229 

0230 

0231 

0232 

0233 

0234 

0235 

0236 

0237 

0238 

0239 

0240 

0241 

0242 

0243 

0244 

0245 

0246 

0247 

0248 

0249 

0250 

0251 

0252 

0253 

0234 

0235 

0256 

0257 

0258 

0239 

0260 

0261 

0262 

0263 

0264 

0265 

0266 

0267 

0268 

0269 

0270 

0271 

0272 

0273 

0274 

0275 

0276 

0277 

0278 

0279 

0280 

0281 

0282 

0283 
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c 

LST«NH-1 

RNK=4.  0*RNH*RNH 

SNK=4.  0*<LST*LST) 

BAC«=BACV 

DAC=RNS 

BAQ=BACV+4  0*RNH*<RNH-1  0) 

C 

C  COMPUTE  ADDRESSES  OF  FOLLOWING  VECTORS 

C  AR=ADDRESS  OF  FIRST  COLUMN  OF  COVARIANCE  MATRIX 

C  AG-ADDRESS  OF  0  COLUMN  OF  COVARIANCE  MATRIX 

C  AC=ADDRESS  OF  FIRST  COLUMN  OF  EXTRACTED  MATRIX 

C 

AR==BAC 

AG=BAG 

AC«DAC 

D  WR ITE (6»  330 )  BAC, BAG,  DAC 

D330  FORMAT ( 1H  , 3F10. 1 ) 

C 

INIT=1 

KK=1 

K=1 

C 

C  CONFIGURE  STEERING  VECTORS  AS  A  SINGLE  COMPLEX  BUFFER 

C 

CALL  MPCLB (IBUS2+7,  BASV#  0.  5*BSSV, 1,  1,  0) 

C 

C  CONFIGURE  EACH  ROW  OF  COVARIANCE  MATRIX  ON  BUS  3<  R  OR.  INV(R)) 

C  AS  A  COMPLEX  VECTOR  USING  BUFFER  NUMBERS  FROM  32  ONWARDS 

C 

DO  426  J-1,LST 

CALL  MPCLB  < IBUS3+31+J# AC+4*< J-l ) ,  RNH-1.  0,  1, LST# 0) 

426  CONTINUE 

CALL  BFPRM(LQBF*  NRCHi  XF,  IER ) 

IFCIER  .  NE.  0)  STOP  1 
C 

C  PREPARE  TO  WRITE  FIRST  BLOCK  ON  MAG.  TAPE 

C 

777  I HEAD  C 1 )  =  ! 

IB1<1 >®NH 
IB1<2)=NBMS 
IB1 (3)“NF 
DO  78  J**l  #  NF 
RFQC J)=SF< J) 

78  CONTINUE 

RFQ<  NF+1 )«D 
RFG(NF+2)=1500.  0 
C 

C  SWITCH  M.  T.  LU#  FROM  7  TO  10  OR  VICE  VERSA 

C 

LUMT  « 1 7— LUMT 

CALL  EXEC (3,  LUMT+400B) 

CALL  EXEC (2# LUMT ,  IHEADi  527) 

I HEAD ( 1 ) =2 
C 
C 

C  CLEAR  THE  COVARIANCE  MATRIX  BEFORE  THE 

C  NEXT  CYCLE  OF  ACCUMULATION  AND  PROCESSING 

C 

77  CALL  MPCLB  ( 20,  BACV,  BSC V,  0,  1,  0> 

CALL  VCLR ( 20 ) 

C 

C  TT=0.  0 

DO  80  1=1#  NA 

C  CALL  EXEC ( Hi  ITIME) 

C  CALL  MBC0V(7, INITi XF#  I END ) 

C  IF ( IEND  .  LT.  0)  GOTO  99 

CALL  MBC2(20i  LOBF#  NRCH) 

C  CALL  EXEC<11, ITIME(6) ) 

C  I D I FF“ I T I ME  <  8 )  —  I T I ME ( 3 ) 

C  IF< IDIFF  .  LT.  0)  IDIFF=IDIFF+60 

C  I D  IFF**  I D I  FF*60+ 1 T I  ME  ( 7 )  —  I T I  ME  <  2 ) 

C  RDIFFt=IDIFF*+ ( I  TIME  (6  >-ITIME<  1 )  >*#0.  01 

C  TT»TT+RDIFF 
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0284 

0285 

0286 

0287  80 

0288  C 
0289  C81 

0290  C 
0291  C 
0292  C 
0293  C 
0294  C 
0295  C 
0296 
0297 
0298 
0299 
0300  C 
0301  C 
0302  C 
0303  C 
0304 
0305 
0306  C 
0307  C 
0308  C 
0309 
0310 
0311 
0312 
0313 
0314 
0315  C 
0316  85 

0317  C 
0318  C 
0319  0 

0320 
0321  C 
0322  C 
0323  C 
0324  C 
0325  C 
0326  C 
0327 
0328 
0329 
0330 
0331 
0332  C 
0333  C 
0334  C 
0335  C 
0336  C 
0337  C 
0338  C 
0339 
0340 
0341 
0342  C 
0343  C 
0344  C 
0345 
0346  C 
0347  D 
0348  D41 5 

0349 
0350  C 
0351  C 
0352  C 
0353 
0354  C 
0355 
0356  C 


CALL  C0VAR(NS# NH, K, NDXSF, NF# XF) 

K*0 
INIT»0 
CONTINUE 
WRITE<6.  81)  TT 

FORMAT (  '  M.  T.  READ  TIME  =  '#  F8.  2) 

CALL  EXEC (11/ IT I ME ) 

CONFIGURE  COVARIANCE  MATRIX  AS  A  REAL  VECTOR  THEN 
NORMALIZE  &  CONJUGATE  BY  NEGATING  THE  ODD  PART 

CALL  MPCLB (20,  BACV, BSCV#  0#  1#  0) 

CALL  VSMAK20,  54,  20,  0) 

CALL  MP0DD(8,  20) 

CALL  VNEG< 8 ) 

CONFIGURE  FIRST  COLUMN  OF  COVARIANCE  MATRIX  AS  A  REAL  VECTOR 
SO  COLUMNS  CAN  BE  MOVED  ONE  BY  ONE  TO  BUS  3  USING  ,,VMOVM 

CALL  MPCLB  <  I BUS2+8,  BACV,  2.  0*RNH,  0,  1,0) 

CALL  MPCLB < IBUS3+9# RNS#  2  0*LST# 0,  1, 0) 

CREATE  FUNCTION  LIST  TO  EXTRACT  A  COLUMN  OF  A  COVARIANCE  MATRIX 

IF <KK  .  EG.  0)  GOTO  85 
CALL  MPBFL(l) 

CALL  VM0V(9i  8) 

CALL  MPSBA(8#  2*NH) 

CALL  MPSBA<9#  2*LST) 

CALL  MPEFL(l) 

CALL  MPFOR  ( 1 #  1$ LST #  1 ) 

PRINT  THE  COVARIANCE  MATRIX  JUST  TRANSFERRED 
I F ( I SSW ( 0 )  LT  0)  CALL  C VMPR (31# LST #  AC  #  SF ( 1 ) ) 

CONFIGURE  THE  FOLLOWING  BUFFERS 

21  *  RESULT  OF  INV(R)  Q(COMPLEX) 

22  «  RESULT  OF  R.B<COMPLEX)  OR  INV ( R >. B ( COMPLEX ) 

23  =  RESULTS  MATRIX (REAL) 

BA=10. 0*NBMS 
STEP«4.  0*( RNH-1.  0) 

CALL  MPCLB (21 » BA# RNH-1 .  0,  1#  1#  0) 

CALL  MPCS0<22, 58# NH-1, 1) 

CALL  MPCLM(23#  2.  0*NBMS#  4#  NBMS) 

CONFIGURE  LOGICAL  BUFFER  CONSISTING  OF  REAL  PARTS 
OF  THE  DIAGONAL  ELEMENTS< EXCEPT  LAST  ONE) 

OF  THE  COVARIANCE  MATRIX 

THEN  COMPUTE  AVERAGE  SUM  IN  SCALAR  132 

AND  TAKE  RECIPROCAL  OF  THIS  VALUE 

CALL  MPCLB (20# BACV,  RNH-1.  0,  0#  2#(NH+1),  0) 

CALL  SSUMU32#  20#  51  ) 

CALL  SDIVU29#  1#  132) 

CONFIGURE  FIRST  STEERING  VECTOR 

CALL  MPCLB  (25,  BASV,  RNH-1.  0#  1#  1,  0) 

WRITE (6» 415) 

FORMAT (1 HI) 

CALL  MPCOL (26#  1, 23) 

CONFIGURE  Q-VECTOR( COLUMN) 

CALL  MPCLB  (24#  AG,  RNH-1.  0#  1#  1#  0) 

IF ( ITYPE  .  LE  2)  GOTO  425 
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0357 
0358 
0359  C 
0360  C 
0361  C 
0362  C 
0363 
0364 
0365  C 
0366  C 
0367  C 
0368  C 
0369 
0370 

0371  420 

0372 
0373 
0374  C 
0375  C 
0376  C 
0377  C 
0378 
0379 
0380 
0381  C 
0382  C 
0383  C 
0384  425 

0385 
0386 
0387 
0388  C 
0389  C 
0390  C 
0391 
0392  C 
0393 
0394 
0395  C 
0396  C 
0397  C 
0398 
0399  C 
0400 
0401  C 
0402  C 
0403  C 
0404  C 
0405 
0406 

0407  430 

0408  C 
0409  C 
0410  C 
0411 
0412  C 
0413  C 
0414  C 
0415  C 
0416  C 
0417 
0418 
0419 
0420 
0421 
0422 
0423 
0424 
0425  C 
0426  435 

0427  C 
0428  C 


MAP*“MPCLA(  IBUS3+31,  AC,  LST,  LST,  1,  0) 

CALL  CMINV(31 # 124, 2) 

CMINV  GIVES  -INV(R)  SO  CONFIGURE  THE 
AREA  AS  A  REAL  BUFFER  Sc  MULTIPLY  BY  -1 

CALL  MPCLB ( 31 , AC, 2  0*(LST*LST>, 0, 1,0) 

CALL  VSMA1 (31,  3,  31, 0) 

COMPUTE  INV (R ) .  Q  IN  SCALAR  TABLE 
AND  SAVE  RESULT  IN  BUFFER  NO.  21 

DO  420  J-1,LST 

CALL  CSD0T(54+2*J,  J+31,  52, 24) 

CONTINUE 

CALL  MPEQBC27,  21) 

CALL  MPTSB ( 27,  56, LST) 

CONFIGURE  ROW  VECTOR  Q*  AND 
COMPUTE  G*  INV<R ) .  G  IN  SCALARS  56-57 

BA^AR+STEP 

CALL  MPCLB<27,  BA,  RNH-1.  0,  1,  NH,  0) 

CALL  CSD0T<56, 27, 52,21) 

EXTRACT  POWER  OF  REFERENCE  PHONE  TO  SCALAR  126 
CONTINUE 

CALL  MPEGB ( 27 ,  24 ) 

CALL  MPSBA(27,  LST) 

CALL  MPTBS ( 27 ,  126,  1 ) 

COMPUTE  G-G*. INV(R). Q  IN  ^CALAR  127 

CALL  SSUB< 127,  126,  56) 

IFCKK  .EG.  0)  GOTO  448 
KK=0 


CREATE  FUNCTION  LIST  TO  PERFORM  PROCESSING  FOR  ONE  BEAM 

CALL  MPBFL< 2 > 

IF ( ITYPE  . GE.  3)  GOTO  435 

CONFIGURE  FIRST  ROW  OF  COVARIANCE  MATRIX 
COMPUTE  R. B  IN  SCALAR  TABLE 

DO  430  J**l,  LST 

CALL  CSDOT< 56+2* J, J+31,  52, 25) 

CONTINUE 

TAKE  CONJUGATE  OF  STEERING  VECTOR 
CALL  CVCNJ ( 25 ) 

COMPUTE  B*.  R  B  IN  SCALARS  120-121 

COMPUTE  B*.  G  IN  SCALARS  56-57 

THEN  ( SCLR  <  56 ) **2+SCLR ( 57 ) **2 ) /SCLR (126) 

CALL  CSDOT( 120,  25,  52, 22) 

CALL  CSD0T(56, 25,  52, 24) 

CALL  SMUL(56,  56,  56) 

CALL  SMUL( 57,  57,  57) 

CALL  SADD(56,  56,  57) 

CALL  SDI V( 56,  56,  126) 

CALL  SSUB ( 121 ,  120,  56) 

GOTO  445 

CONTINUE 

COMPUTE  INV(R ) .  B  IN  SCALAR  TABLE 
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<- 


0429 

0430 

0431 

0432 

0433 

0434 

0435 

0436 

0437 

0438 

0439 

0440 

0441 

0442 

0443 

0444 

0445 

0446 

0447 

0448 

0449 

0450 

0451 

0452 

0453 

0454 

0455 

0456 

0457 

0458 

0459 

0460 

0461 

0462 

0463 

0464 

0465 

0466 

0467 

0468 

0469 

0470 

0471 

0472 

0473 

0474 

0475 

0476 

0477 

0478 

0479 

0480 

0481 

0482 

0483 

0484 

0485 

0486 

0487 

0488 

0489 

0490 

0491 

0492 

0493 

0494 

0495 

0496 

0497 

0498 

0499 

0500 


C 


440 

C 

C 

C 


C 


C 

445 

C 

C 

c 

c 

c 


c 

448 


C 

C 

C 

C 

C 


DO  440  J-l.LST 

CALL  CSD0T(56+2*J,  0+31,  52,  25) 

CONTINUE 

TAKE  CONJUGATE  OF  STEERING  VECTOR 
CALL  CVCNJ ( 25 ) 

COMPUTE  B*.  INV(R)  B  IN  SCALARS  122-123 
COMPUT  B*  INV(R).Q  IN  SCALARS  56~57 
THEN  ( SCLR ( 56  >  **2+SCLR ( 57 ) **2 ) /SCLR (127) 

CALL  CSDOT ( 122#  25. 52,  22) 

CALL  CSDOT ( 56*  25,  52,  21 ) 

CALL  SMUL<56.  56,  56) 

CALL  SMUL<57,  57.  57) 

CALL  SADDC56,  56.  57) 

CALL  SDI V ( 56.  56,  127) 

CALL  SADD( 123,  122,  56) 

CALL  SDIVU22,  50,  122) 

CALL  SDIV< 123,  50,  123) 

CONTINUE 

OUR  4  RESULTS  ARE  NOW  IN  SCALARS  120-123 
SO  TRANSFER  THEM  TO  NEXT  COLUMN  OF  MATRIX 

CALL  MPTSB<26, 120, 4) 

CALL  MPSBA(25, LST) 

CALL  MPEFL < 2 ) 

CALL  MPIBC(O) 

CALL  MPFOR  <1,1,  NBMS,  2 ) 

CALL  MPRBC(O) 

AFTER  EXECUTION  OF  FUNCTION  LIST  2  FOR  EACH  BEAM 
ALL  THE  STEERING  VECTORS  HAVE  BEEN  CONJUGATED 
SO  RESTORE  THEM  READY  FOR  THE  NEXT  CYCLE 


CALL  CVCNJ <  7 ) 

C 

C  PROCESSING  FOR  SELECTED  FREQUENCY  COMPLETE 

C  FOR  SELECTED  COLUMN  OF  RESULTS  MATRIX 

C  MULTIPLY  RESULTS  BY  1.0/<AVERAGE  POWER) 

C  FIND  MAXIMUM  VALUE  AND  COMPUTE  1.0/MAX 

C  THEN  MULTIPLY  EACH  ELEMENT  BY  THIS  VALUE 

C  CONVERT  THESE  VALUES  TO  DECIBELS, 

C  TRANSFER  TO  HOST  &  PLOT 

C 

CALL  MPR0W<26, I TYPE, 23) 

CALL  VSMA1 (26,  129, 26, 0) 

CALL  SMAX ( 131 ,  26,  127) 

CALL  MPRST < 131 , PMAX, 2,  1 ) 

CALL  SDIV( 127, 1, 131 ) 

CALL  VSMAK26,  127,  26,  0) 

CALL  VL0GH(26, 31, 26) 

CALL  MPRDB <26, POWR (4) ,  1,  1 ) 

J-3+NBMS 
DO  449  1-4,  J 

POWR  < I ) -RNORM (POWR ( I ) )-5.  0 
449  CONTINUE 

CALL  EXEC ( 20, 0, POWR, ILEN, 0, 0, ICLS) 

C  CALL  EXEC (lit  I T I ME < 6 ) ) 

C  WRITE(6, 446)  < ITIME < J ) . J=10, 1, -1 ) 

C446  FORMAT <  '  TIME:  ',1016) 

C 


IF < IFBRK(O) )  90,77 
90  WRITE( 1,  452) 

452  FORMAT ( '  EX  *  EXIT  PROGRAM  ' / 9  ST  =  SWITCH  TAPES') 
READ( 1,453)  IANS 
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0501 

453 

FORMAT (A2) 

0502 

IF( IANS  .  EQ.  IEX )  GOTO  99 

0503 

IF (IANS  .  NE.  1ST)  GOTO  90 

0504 

C 

0305 

C 

WRITE  EOF  ON  M.  T.  &  REWIND 

0506 

c 

0507 

CALL  EXEC ( 3, LUMT + 1 OOB ) 

0508 

CALL  EXEC <3, LUMT+400B) 

0509 

GOTO  777 

0510 

99 

CALL  MPCLS(O) 

0311 

CALL  EXEC<3,  LUMT+10QB) 

0512 

NPTS«-1 

0513 

CALL  EXEC(20,  0,  POWR,  2,  0,  0, 

0514 

CALL  EXEC (6) 

0515 

END 

0516 

END* 

* 


&CVM8I  T=0OOO4  IS  ON  CR00013  USING  00017  BLKS  R=0000 


0001 

0002 

0003 

0004 

0003 

0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 

0022 

0023 

0024 

0023 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 


FTN4,  L 

SUBROUTINE  COVAR(NS, NH, K, ISF, NSF, XF) 
DIMENSION  ISF (50) 

COMMON  LUMT/ IHEAD ( 1 5 > * HRB (800 ) 

DATA  IBUS2/128/,  IBUS3/192/ 

IF ( K  .  EG.  0)  GOTO  100 
RNS=*NS 
RNH=NH 
BS*RNS*RNH 
CVSIZ»2.  0*'RNH*RNH 
ISIZ=2*NH*NSF 
BACVaa8191 8.  0—2.  0*CVSIZ 
ACM-RNS+4.  0*RNH*NSF 
C 

C  SET  SCALE  FACTOR (XF)  IN  SCALAR  128 

C  RAW  DATA  IS  MULTIPLIED  BY  XF  BEFORE  FFT 

C 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


C 

C 

C 

C 


CALL  MPWST ( 128 1  XF,  1/  1 ) 

CONFIGURE  THE  FOLLOWING  PERMANENT  BUFFERS 

1  *  CONSTANT  COMPLEX  VECTOR ( CURRENT  VALUE  OF  SCALARS  56-57) 

2  -  DATA  AREA  AS  A  SINGLE  REAL  BUFFER 

29  =  DATA  AREA  AS  A  REAL  MATRIX 

10  «  SOURCE  BUFFER (REAL)  FOR  FFT 

11  *  COSINE  TABLE  FOR  FFT 

12  »  COMPLEX  BUFFER  FOR  RESULT  OF  FFT 

30  =  MATRIX  OF  SELECTED  FREQUENCIES 

28  ■  COVARIANCE  MATRIX  FOR  SELECTED  FREQUENCY 

16  ■  FIRST  COLUMN  OF  COVARIANCE  INCREMENT  MATRIX 

17  *  COVARIANCE  INCREMENT ( CONFIGURED  AS  REAL  VECTOR) 

CALL  MPDCV ( 1 i 56, NH, 1 ) 

CALL  MPCLB  ( IBUS2+2,  0.  0,  BS,  0,  1, 0) 

CALL  MPCLM(29, 0.  0,  NH,  NS) 

CALL  MPCLB ( IBUS3+10,  RNS,  RNS,  0,  1,0) 

CALL  MPCLB  ( IBUS3+U,  O.  0,  0  5*RNS,  O,  1,0) 

CALL  MPCLB ( IBUS3+12,  RNS,  0.  5*RNS,  1,  1,  0) 

CALL  MPCLM(30, RNS,  2*NH,  NSF) 

CALL  MPCLBC28, BACV,  CVSIZ,  0,  1,  O) 

CALL  MPCLB (IBUS3+16, ACM,  RNH,  1,  1,  0) 

CALL  MPCLB (IBUS3+17, ACM,  CVSIZ,  0,  1,  0) 

CREATE  FUNCTION  LIST  14  TO  MOVE  ONE  ROW  TO  BUS  3, 

PERFORM  FFT  AND  MOVE  RESULT  BACK  TO  ORIGINAL  ROW 

CALL  MPBFLU4) 

CALL  VMOV< 10,3) 

CALL  FFTR3 ( 12,  10,  11 ) 

CALL  VMOV ( 3,  10) 

CALL  MPNXR(3,  29) 

CALL  MPEFL (14) 


r 
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« 


% 


0052  C 

0053  C  CREATE  FUNCTION  LIST  15  TO  MOVE 

0054  C  2  COLUMNS  CORRESPONDING  TO  ONE  FREQUENCY 

0055  C 

0056  CALL  MPBFLU5) 

0057  CALL  MPEVN< 15, 13) 

0058  CALL  MPODD ( 14#  13) 

0059  CALL  VMOV( 15/4) 


0060 
0061 
0062 
0063  C 
0064  C 
0065  C 
0066  C 
0067 
0068 
0069 
0070 
0071 
0072 

0073  25 

0074 
0075  C 
0076  C 
0077  C 
0078  C 
0079 
0080 
0081 
0082 
0083 
0084  C 
0085  C 
0086  C 
0087  C 
0088 
0089  C 
0090  C 
0091  C 
0092  C 
0093 
0094 
0095 
0096  C 
0097 
0098  C 
0099  C 
0100  C 
0101 
0102  C 
0103  C 
0104  C 
0105  C 
0106 
0107  C 
0108  C 
0109  C 
0110 
0111 
0112  C 
0113  C 
0114  C 
0115  C 
0116 
0117 
0118  C 
0119  C 
0120  C 
0121  C 


CALL  VMOV (14/  5) 

CALL  MPNXC ( 13# 30) 

CALL  MPEFL (15) 

CREATE  FUNCTION  LIST  13  TO  SELECT  S<  MOVE 
TWO  COLUMNS  OF  FREQUENCY  DATA 

CALL  MPBFL ( 13 ) 

DO  25  1*1/ NSF 
J=*2* ISF ( I  )  —1 
CALL  MPCOL  (  4/  J,  29) 

CALL  MPCOL<  5/ J+li 29) 

CALL  MPXFL ( 1 5 ) 

CONTINUE 
CALL  MPEFL (13) 

CREATE  FUNCTION  LIST  12  TO  COMPUTE  INCREMENT  FOR 
ONE  COLUMN  OF  A  COVARIANCE  MATRIX 

CALL  MPBFL (12) 

CALL  MPTBS(19, 56, 1 ) 

CALL  CCVML(15,  0,  1/  0,  18) 

CALL  MPSBA ( 1 5/ NH ) 

CALL  MPEFL (12) 

CREATE  FUNCTION  LIST  11  TO  CALCULATE  THE 
INCREMENTS  Sc  UPDATE  A  COVARIANCE  MATRIX 

CALL  MPBFL (11) 

MULTIPLY  SELECTED  COLUMN  BY  COMPLEX  CONJUGATE  OF  EACH  ELEMENT 
AND  SAVE  THE  VECTOR  SO  FORMED  IN  BUFFER  15 

CALL  MPEQB( 15,  16) 

CALL  MPFORd,  1,  NH,  12) 

CALL  VAD(20,  17) 

CALL  MPEFL (11) 

CREATE  FUNCTION  LIST  9  TO  PERFORM  COMPLETE  PROCESSING 
CALL  MPBFL (9) 

RE-DEFINE  THE  FOLLOWING  BUFFER 

18  *  FIRST  COLUMN  OF  FREQUENCY  MATRIX 

CALL  MPCLB ( IBUS3+18,  RNS,  RNH,  1,  1,  O) 

CONVERT  DATA  TO  REAL  AND  SCALE 

CALL  VFLTC2,  39,  20,  O) 

CALL  VSMA1 (2,  128, 2, 0) 

DEFINE  FIRST  ROW  OF  DATA 
Sc  PERFORM  FFT  ON  EACH  ROW 

CALL  MPROW ( 3,  1, 29) 

CALL  MPFORCl,  1 , NH,  14) 

DEFINE  FIRST  COLUMN  OF  SELECTED  FREQUENCIES 
AND  EXTRACT  THE  SELECTED  FREQUENCY  DATA 
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0122 

CALL  MPCOL ( 13/  1,  30) 

0123 

CALL  MPXFL ( 13 ) 

0124 

C 

0125 

C 

SELECT  FIRST  COVARIANCE  MATRIX 

0126 

C 

AND  UPDATE  EACH  COVARIANCE  MATRIX 

0127 

C 

0128 

CALL  MPEGB(19,  18) 

0129 

CALL  MPXFL (11) 

0130 

CALL  MPEFL ( 9 ) 

0131 

C 

0132 

C 

REDEFINE  DATA  AREA  AS  AN  INTEGER  BUFFER 

0133 

C 

0134 

100 

CALL  MPCLB  ( IBUS2+20,  0.  0,  BS,  2,  2,  0> 

0135 

CALL  MPXFL (9) 

0136 

CALL  MPCLB (31  i  RNS,  1. 0*ISIZ, 0.  1* 0) 

0137 

CALL  MPRDBOl,  HRB,  1,1) 

0138 

DO  150  J=l,  ISIZ 

0139 

HRB  <  J  >  *RN0RM ( HRB  <  J ) ) 

0140 

150 

CONTINUE 

0141 

CALL  EXEC (2, LUMT , I  HEAD, 2*ISIZ+15) 

0142 

I HEAD  < 1 >  =  I HEAD <  1 )  + 1 

0143 

C 

0144 

RETURN 

0145 

END 

0146 

END* 
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